A small OSL libray

The code below is meant to be saved as `equations.h` in a location that is looked in by OSLs `#include` directive. Currently is contains the following two functions (but more will be added as I develop other shaders):

• `cubic(float A, float X, int L)`
• `point cubicspline(float t, point P0, point P1, point P2, point P3)`
• `int splinedist(point p0, point p1, point p2, point Pos, float d, float tc)`
The first one will find the real roots of the cubic equation
```A * t**3 + A * t**2 + A * t + A == 0
```
and will return the roots in `X` and the number of (real) roots in `L`. I converted the code from http://van-der-waals.pc.uni-koeln.de/quartic/quartic.html to OSL, which is supprisingly simple. There is other cool stuff on that site, so check it out.

The second one will return a point on a cubic bezier spline that starts at P0, heads in the direction of P1 and ends at P3 coming from the direction of P2. (note that the `spline()` only works for uniformly distributed points and the documentation is a bit unclear so I thought I'd better implement one I understood)

The third one will return the closest distance to a cubic spline with control points P0, P1 and P2, returning the result in d and the corresponding t value in tc. The function returns 1 if a closest distance is found with a t value in the range [0,1], zero otherwise.

An example of how to include this code is shown below

```#include "equation.h"

{
...
point p = cubicspline(t,P0,P1,P2,P3);
...
}
```

equations.h

This is the actual library. Save it for example as `scripts\addons\cycles\shaders\equations.h` in your Blender installation directory or in the same directory as your main .osl file.

```// cubic roots adapted for OSL from http://van-der-waals.pc.uni-koeln.de/quartic/quartic.html
// so now we have a translation from Fortran -> C -> OSL. It doesn't look that well, but it works :-)

float CBRT(float Z) { return abs(pow(abs(Z),1.0/3.0)) * sign(Z); }

/*-------------------- Global Function Description Block ----------------------
*
*     ***CUBIC************************************************08.11.1986
*     Solution of a cubic equation
*     Equations of lesser degree are solved by the appropriate formulas.
*     The solutions are arranged in ascending order.
*     ******************************************************************
*     A(0:3)      (i)  vector containing the polynomial coefficients
*     X(1:L)      (o)  results
*     L           (o)  number of valid solutions (beginning with X(1))
*     ==================================================================
*   17-Oct-2004 / Raoul Rausch
*  Conversion from Fortran to C
*
*-----------------------------------------------------------------------------
*/

int cubic(float A, float X, int L)
{
float PI = 3.1415926535897932;
float THIRD = 1./3.;
float U,W, P, Q, DIS, PHI;
int i;

// ====determine the degree of the polynomial ====

if (A != 0.0)
{
//cubic problem
W = A/A*THIRD;
P = pow((A/A*THIRD - pow(W,2)),3);
Q = -.5*(2.0*pow(W,3)-(A*W-A)/A );
DIS = pow(Q,2)+P;
if ( DIS < 0.0 )
{
//three real solutions!
//Confine the argument of ACOS to the interval [-1;1]!
PHI = acos(min(1.0,max(-1.0,Q/sqrt(-P))));
P=2.0*pow((-P),(5.e-1*THIRD));
for (i=0;i<3;i++)
U[i] = P*cos((PHI+2*((float)i)*PI)*THIRD)-W;
X = min(U, min(U, U));
X = max(min(U, U),max( min(U, U), min(U, U)));
X = max(U, max(U, U));
L = 3;
}
else
{
// only one real solution!
DIS = sqrt(DIS);
X = CBRT(Q+DIS)+CBRT(Q-DIS)-W;
L=1;
}
}
else if (A != 0.0)
{
P = 0.5*A/A;
DIS = pow(P,2)-A/A;
if (DIS > 0.0)
{
// 2 real solutions
X = -P - sqrt(DIS);
X = -P + sqrt(DIS);
L=2;
}
else
{
// no real solution
L=0;
}
}
else if (A != 0.0)
{
//linear equation
X =A/A;
L=1;
}
else
{
//no equation
L=0;
}
/*
*     ==== perform one step of a newton iteration in order to minimize
*          round-off errors ====
*/
for (i=0;i < L;i++)
{
X[i] = X[i] - (A+X[i]*(A+X[i]*(A+X[i]*A)))
/(A+X[i]*(2.0*A+X[i]*3.0*A));
}

return 0;
}

point cubicspline(float t, point P0, point P1, point P2, point P3)
{
return (1-t)*(1-t)*(1-t)*P0
+ 3*(1-t)*(1-t)*t*P1
+ 3*(1-t)*t*t*P2
+ t*t*t*P3;

// how to group all factors of the powers of t
//
// (1-2t+t2-t+2t2+t3)P0 => (1-3t+3t2+t3)P0
// (3t-6t2+3t3)P1
// (3t2-3t3)P2
// t3P3
//
// 1 * ( P0 )
// t * ( 3P0+3P1 )
// t2* ( 3P0-6P1+3P2 )
// t3* ( P0+3P1-3P2+P3)
}

// calculate the closest distance from Pos to a quadratic bezier curve
// the curve is defined by the 3 points P0, P1 and P2

int splinedist(point p0, point p1, point p2, point Pos, float d, float tc){

point P0 = p0;
point P1 = p1;
point P2 = p2;

// following definitions are for the four polynomic coefficients for the well known
// equation dB/dt . (Pos-B)  (i.e. the inproduct of the tangent to the bezier and the
// difference vector from the point under considertion to the Bezier curve.
// If the difference vector is perpendicular to the tangent we have found a closest point
// on the Bezier curve.
// The stuff below is generated by a script and no effort is spent on collecting factors.
// We let the OSL compiler worry about that :-)

float t0 = -2*P0*P0+-2*P1*Pos+2*P0*P1+2*P0*Pos
+-2*P0*P0+-2*P1*Pos+2*P0*P1+2*P0*Pos
+-2*P0*P0+-2*P1*Pos+2*P0*P1+2*P0*Pos;

float t1 = -4*P0*P1+-4*P0*P1+-4*P0*P1+-2*P0*Pos+-2*P2*Pos+2*P0*P0+2*P0*P2+4*P0*P0+4*P1*P1+4*P1*Pos
+-4*P0*P1+-4*P0*P1+-4*P0*P1+-2*P0*Pos+-2*P2*Pos+2*P0*P0+2*P0*P2+4*P0*P0+4*P1*P1+4*P1*Pos
+-4*P0*P1+-4*P0*P1+-4*P0*P1+-2*P0*Pos+-2*P2*Pos+2*P0*P0+2*P0*P2+4*P0*P0+4*P1*P1+4*P1*Pos;

float t2 = -8*P1*P1+-4*P0*P0+-4*P0*P2+-4*P1*P1+-2*P0*P0+-2*P0*P2+2*P0*P1+2*P1*P2+4*P0*P1+4*P0*P1+4*P1*P2+8*P0*P1
+-8*P1*P1+-4*P0*P0+-4*P0*P2+-4*P1*P1+-2*P0*P0+-2*P0*P2+2*P0*P1+2*P1*P2+4*P0*P1+4*P0*P1+4*P1*P2+8*P0*P1
+-8*P1*P1+-4*P0*P0+-4*P0*P2+-4*P1*P1+-2*P0*P0+-2*P0*P2+2*P0*P1+2*P1*P2+4*P0*P1+4*P0*P1+4*P1*P2+8*P0*P1;

float t3 = -4*P0*P1+-4*P0*P1+-4*P1*P2+-4*P1*P2+2*P0*P0+2*P0*P2+2*P0*P2+2*P2*P2+8*P1*P1
+-4*P0*P1+-4*P0*P1+-4*P1*P2+-4*P1*P2+2*P0*P0+2*P0*P2+2*P0*P2+2*P2*P2+8*P1*P1
+-4*P0*P1+-4*P0*P1+-4*P1*P2+-4*P1*P2+2*P0*P0+2*P0*P2+2*P0*P2+2*P2*P2+8*P1*P1;

float A = {t0,t1,t2,t3};
float T ;
int n ;
cubic(A,T,n);
d = 1e6;
// cubic() will return 0 , 1 or 3 values for t
// we are only interested in values that lie in the interval [0,1]
// for those we calculate the position on the curve and check whether
// we have found the shortest distance.
int found = 0;
while(n>0){
n--;
if(T[n]>=0 && T[n]<=1){
float t = T[n];
found = 1;
float dd = distance((1-t)*(1-t)*P0 + 2*(1-t)*t *P1 + t*t*P2, Pos);
if (dd < d) {
d = dd;
tc = t;
}
}
}
return found;
}
```