#define TRACE_CUBIC

#ifdef TRACE_CUBIC
#include <math.h>
#endif

#ifdef TRACE_CUBIC
/*************************************************************************/
/*                                                                       */
/* Returns the one-dimensional position of the Bezier point with         */
/* parametric value t and control point values s0 .. s3.                 */
/*                                                                       */
static double bez_param(double s0, double s1, double s2, double s3, double t)
{
   double u = 1 - t;

   return   u * u * u * s0
          + 3 * t * u * (u * s1 + t * s2)
          + t * t * t * s3;
}


/*************************************************************************/
/*                                                                       */
/* Returns the unnormalised s-coordinate of the point p in the r-s       */
/* coordinate system for the Bezier curve defined by the points          */
/* arc[0..3].                                                            */
/*                                                                       */
static FT_Pos s_dist(FT_Vector* arc, FT_Vector* p)
{
   return   (arc[0].y - p->y) * (arc[3].x - arc[0].x) 
          - (arc[0].x - p->x) * (arc[3].y - arc[0].y);
}


/*************************************************************************/
/*                                                                       */
/* Returns the distance between points p and q.                          */
/*                                                                       */
static double dist(FT_Vector* p, FT_Vector* q)
{
   FT_Pos dx = q->x - p->x;
   FT_Pos dy = q->y - p->y;

   return sqrt(dx * dx + dy * dy);
}


/*************************************************************************/
/*                                                                       */
/* Returns the exact maximum lateral deviation of the Bezier curve       */
/* defined by the points arc[0..3].                                      */
/*                                                                       */
static double max_dev(FT_Vector* arc)
{
   FT_Pos s1 = s_dist(arc, &arc[1]);
   FT_Pos s2 = s_dist(arc, &arc[2]);

   double v, t_max, s_max;

   if (s1 == 0 && s2 == 0)
      return 0;

   v = abs(s1) < abs(s2) ? (double)s1 / s2 : (double)s2 / s1;

   t_max = v == 1 ? 0.5 : (2 - v - sqrt(v * v - v + 1)) / (3 * (1 - v));

   s_max = bez_param(0, 1, v, 0, t_max);

   return (abs(s1) < abs(s2) ? abs(s2) : abs(s1)) * s_max / dist(&arc[0], &arc[3]);
}
#endif /* TRACE_CUBIC */


  static void
  gray_render_cubic( RAS_ARG_ const FT_Vector*  control1,
                              const FT_Vector*  control2,
                              const FT_Vector*  to )
  {
#ifdef TRACE_CUBIC
    static int    tot_arcs = 0;
    static int    tot_segs = 0;
    static double tot_dev = 0;
#endif

    TPos        dx, dy, da, db;
    int         top, level;
    int*        levels;
    FT_Vector*  arc;


    dx = DOWNSCALE( ras.x ) + to->x - ( control1->x << 1 );
    if ( dx < 0 )
      dx = -dx;
    dy = DOWNSCALE( ras.y ) + to->y - ( control1->y << 1 );
    if ( dy < 0 )
      dy = -dy;
    if ( dx < dy )
      dx = dy;
    da = dx;

    dx = DOWNSCALE( ras.x ) + to->x - 3 * ( control1->x + control2->x );
    if ( dx < 0 )
      dx = -dx;
    dy = DOWNSCALE( ras.y ) + to->y - 3 * ( control1->x + control2->y );
    if ( dy < 0 )
      dy = -dy;
    if ( dx < dy )
      dx = dy;
    db = dx;

    level = 1;
    da    = da / ras.cubic_level;
    db    = db / ras.conic_level;
    while ( da > 0 || db > 0 )
    {
      da >>= 2;
      db >>= 3;
      level++;
    }

#ifdef TRACE_CUBIC
    arc      = ras.bez_stack;
    arc[0].x = UPSCALE( to->x );
    arc[0].y = UPSCALE( to->y );
    arc[1].x = UPSCALE( control2->x );
    arc[1].y = UPSCALE( control2->y );
    arc[2].x = UPSCALE( control1->x );
    arc[2].y = UPSCALE( control1->y );
    arc[3].x = ras.x;
    arc[3].y = ras.y;
#endif

    if ( level <= 1 )
    {
      TPos   to_x, to_y, mid_x, mid_y;

      to_x  = UPSCALE( to->x );
      to_y  = UPSCALE( to->y );
      mid_x = ( ras.x + to_x +
                3 * UPSCALE( control1->x + control2->x ) ) / 8;
      mid_y = ( ras.y + to_y +
                3 * UPSCALE( control1->y + control2->y ) ) / 8;

#ifdef TRACE_CUBIC
      gray_split_cubic( arc );
      tot_dev += max_dev(arc);
      tot_dev += max_dev(arc + 3);
      tot_segs+=2;
#endif

      gray_render_line( RAS_VAR_ mid_x, mid_y );
      gray_render_line( RAS_VAR_ to_x, to_y );
      goto Cubic_Done;
    }

#ifndef TRACE_CUBIC
    arc      = ras.bez_stack;
    arc[0].x = UPSCALE( to->x );
    arc[0].y = UPSCALE( to->y );
    arc[1].x = UPSCALE( control2->x );
    arc[1].y = UPSCALE( control2->y );
    arc[2].x = UPSCALE( control1->x );
    arc[2].y = UPSCALE( control1->y );
    arc[3].x = ras.x;
    arc[3].y = ras.y;
#endif

    levels    = ras.lev_stack;
    top       = 0;
    levels[0] = level;

    while ( top >= 0 )
    {
      level = levels[top];
      if ( level > 1 )
      {
        /* check that the arc crosses the current band */
        TPos  min, max, y;

        min = max = arc[0].y;
        y = arc[1].y;
        if ( y < min ) min = y;
        if ( y > max ) max = y;
        y = arc[2].y;
        if ( y < min ) min = y;
        if ( y > max ) max = y;
        y = arc[3].y;
        if ( y < min ) min = y;
        if ( y > max ) max = y;
        if ( TRUNC( min ) >= ras.max_ey || TRUNC( max ) < 0 )
          goto Draw;

        gray_split_cubic( arc );
        arc += 3;
        top ++;
        levels[top] = levels[top - 1] = level - 1;
        continue;
      }

    Draw:
      {
        TPos  to_x, to_y, mid_x, mid_y;

        to_x  = arc[0].x;
        to_y  = arc[0].y;
        mid_x = ( ras.x + to_x + 3 * ( arc[1].x + arc[2].x ) ) / 8;
        mid_y = ( ras.y + to_y + 3 * ( arc[1].y + arc[2].y ) ) / 8;

#ifdef TRACE_CUBIC
        gray_split_cubic( arc );
        tot_dev += max_dev(arc);
        tot_dev += max_dev(arc + 3);
        tot_segs += 2;
#endif

        gray_render_line( RAS_VAR_ mid_x, mid_y );
        gray_render_line( RAS_VAR_ to_x, to_y );
        top --;
        arc -= 3;
      }
    }

Cubic_Done:

#ifdef TRACE_CUBIC
    tot_arcs++;
    printf("%d source arcs: %d segs, avg segs/arc = %f, avg dev/seg = %f\n", 
               tot_arcs, tot_segs, (float)tot_segs / tot_arcs, tot_dev / tot_segs);
#endif

    return;
  }