/* 
 *  Given Bezier curve:
 *
 *    |               |
 *    A          C    |
 *    |*.      .''**-.D
 *    |  *----*       |
 *    0---------B-----1
 *
 * With the constraints:
 *    Ax = 0
 *    Dx = 1
 *    0 <= Bx <= 1
 *    0 <= Cx <= 1
 *
 * It has the parametric formula for point P along the curve, at time 0<=t<=1:
 *    P(t) = t^3 (1+3B-3C) + t^2(-6B+3C) + t(3B)
 *
 * We would like to have an explicit function: Py(x)
 *
 * This requires rearranging, solving the cubic, then passing it to the 
 * parametric function. 
 *
 *    F(t,x) = t^3(1+3Bx-3Cx) + t^2(-6Bx+3Cx) + t(3Bx) - x
 *    solve F(t,x) = 0 for t.
 *
 * A full cubic rootfind will waste time deciding on if there is 1 or 3 real 
 * roots, involving trig or complex numbers etc.
 *
 * We're only looking for the one (guaranteed) root within the unit interval.
 * I've chosen to use the Illinois method which is a bracketed option, and 
 * makes sure we get a real value without risk of diverging.
 *
 *    t = Illinois( F, 0, 1 )
 *
 * Since our bracket starts as the unit interval, the first iteration is very
 * quick to compute the cubics at each end, since t0 is 0, and t1 is 1. Few
 * iterations are needed, no more than 10 for our purposes of animation curves,
 * for sane curves, it will early exit around 3-4.
 *
 */

f32 explicit_bezier( f32 A[2], f32 B[2], f32 C[2], f32 D[2], f32 x )
{
   f32 dAxDx = D[0]-A[0],
       unitBx = (B[0] - A[0]) / dAxDx,
       unitCx = (C[0] - A[0]) / dAxDx,

       /* cubic coefficients */
       a =  3.0f*unitBx - 3.0f*unitCx + 1.0f,
       b = -6.0f*unitBx + 3.0f*unitCx,
       c =  3.0f*unitBx,
       d = -(x - A[0]) / dAxDx,

        t0 = 0.0f,
       Ft0 = d,
        t1 = 1.0f,
       Ft1 = a+b+c+d,
        tc, Ftcx;

   /* Illinois method to find root */
   for( u32 j=0; j<8; j ++ )
   {
      tc = t1 - Ft1*(t1-t0)/(Ft1-Ft0);
      Ftcx = tc*tc*tc*a + tc*tc*b + tc*c + d;

      if( fabsf(Ftcx) < 0.00001f )
         break;
      
      if( Ft1*Ftcx < 0.0f )
      {
         t0 = t1;
         Ft0 = Ft1;
      }
      else
         Ft0 *= 0.5f;

      t1 = tc;
      Ft1 = Ftcx;
   }

   /* Evaluate parametric bezier */
   f32 t2 = tc*tc,
       t3 = tc*tc*tc;
   
   return  D[1] * t3
         + C[1] * (-3.0f*t3 + 3.0f*t2)
         + B[1] * ( 3.0f*t3 - 6.0f*t2 + 3.0f*tc)
         + A[1] * (-1.0f*t3 + 3.0f*t2 - 3.0f*tc + 1.0f);
}