home *** CD-ROM | disk | FTP | other *** search
/ The World of Computer Software / World_Of_Computer_Software-02-387-Vol-3of3.iso / g / gs252src.zip / GS252 / GXPCOPY.C < prev    next >
C/C++ Source or Header  |  1992-09-18  |  17KB  |  502 lines
  1. /* Copyright (C) 1992 Aladdin Enterprises.  All rights reserved.
  2.    Distributed by Free Software Foundation, Inc.
  3.  
  4. This file is part of Ghostscript.
  5.  
  6. Ghostscript is distributed in the hope that it will be useful, but
  7. WITHOUT ANY WARRANTY.  No author or distributor accepts responsibility
  8. to anyone for the consequences of using it or for whether it serves any
  9. particular purpose or works at all, unless he says so in writing.  Refer
  10. to the Ghostscript General Public License for full details.
  11.  
  12. Everyone is granted permission to copy, modify and redistribute
  13. Ghostscript, but only under the conditions described in the Ghostscript
  14. General Public License.  A copy of this license is supposed to have been
  15. given to you along with Ghostscript so you can know your rights and
  16. responsibilities.  It should be in a file named COPYING.  Among other
  17. things, the copyright notice and this notice must be preserved on all
  18. copies.  */
  19.  
  20. /* gxpcopy.c */
  21. /* Path copying and flattening for Ghostscript library */
  22. #include "math_.h"
  23. #include "gx.h"
  24. #include "gserrors.h"
  25. #include "gxfixed.h"
  26. #include "gxarith.h"
  27. #include "gzpath.h"
  28.  
  29. /* Forward declarations */
  30. private int copy_path(P4(const gx_path *, gx_path *, fixed, int));
  31. private int flatten_recur(P8(gx_path *,
  32.   fixed, fixed, fixed, fixed, fixed, fixed, fixed));
  33. private int flatten_sample(P8(gx_path *, int,
  34.   fixed, fixed, fixed, fixed, fixed, fixed));
  35.  
  36. /* Copy a path */
  37. int
  38. gx_path_copy(const gx_path *ppath_old, gx_path *ppath, int init)
  39. {    return copy_path(ppath_old, ppath, max_fixed, init);
  40. }
  41.  
  42. /* Flatten a path. */
  43. /* If flatness is zero, use sampling rather than subdivision: */
  44. /* this is important for Type 1 fonts. */
  45. private fixed
  46. scale_flatness(floatp flatness)
  47. {    fixed scaled_flat = float2fixed(flatness);
  48.     return (scaled_flat > int2fixed(100) ? int2fixed(100) :
  49.         scaled_flat < fixed_half ? fixed_0 :
  50.         scaled_flat);
  51. }
  52. int
  53. gx_path_flatten(const gx_path *ppath_old, gx_path *ppath, floatp flatness, int in_BuildChar)
  54. {    return copy_path(ppath_old, ppath,
  55.              (in_BuildChar ? fixed_0 : scale_flatness(flatness)),
  56.              1);
  57. }
  58.  
  59. /* Add a flattened curve to a path. */
  60. int
  61. gx_path_add_flattened_curve(gx_path *ppath,
  62.   fixed x1, fixed y1, fixed x2, fixed y2, fixed x3, fixed y3,
  63.   floatp flatness)
  64. {    return flatten_recur(ppath, x1, y1, x2, y2, x3, y3,
  65.                  scale_flatness(flatness));
  66. }
  67.  
  68. /* Copy a path, optionally flattening it. */
  69. /* If the copy fails, free the new path. */
  70. private int
  71. copy_path(const gx_path *ppath_old, gx_path *ppath, fixed scaled_flat,
  72.   int init)
  73. {    gx_path old;
  74.     const segment *pseg;
  75.     int code;
  76. #ifdef DEBUG
  77. if ( gs_debug['p'] )
  78.     gx_dump_path(ppath_old, "before copy_path");
  79. #endif
  80.     old = *ppath_old;
  81.     if ( init )
  82.         gx_path_init(ppath, &ppath_old->memory_procs);
  83.     pseg = (const segment *)(old.first_subpath);
  84.     while ( pseg )
  85.        {    switch ( pseg->type )
  86.            {
  87.         case s_start:
  88.             code = gx_path_add_point(ppath, pseg->pt.x, pseg->pt.y);
  89.             break;
  90.         case s_curve:
  91.            {    curve_segment *pc = (curve_segment *)pseg;
  92.             if ( scaled_flat == max_fixed )    /* don't flatten */
  93.                 code = gx_path_add_curve(ppath,
  94.                     pc->p1.x, pc->p1.y,
  95.                     pc->p2.x, pc->p2.y,
  96.                     pc->pt.x, pc->pt.y);
  97.             else
  98.                 code = flatten_recur(ppath,
  99.                     pc->p1.x, pc->p1.y,
  100.                     pc->p2.x, pc->p2.y,
  101.                     pc->pt.x, pc->pt.y,
  102.                     scaled_flat);
  103.             break;
  104.            }
  105.         case s_line:
  106.             code = gx_path_add_line(ppath, pseg->pt.x, pseg->pt.y);
  107.             break;
  108.         case s_line_close:
  109.             code = gx_path_close_subpath(ppath);
  110.             break;
  111.            }
  112.         if ( code )
  113.            {    gx_path_release(ppath);
  114.             if ( ppath == ppath_old )
  115.                 *ppath = old;
  116.             return code;
  117.            }
  118.         pseg = pseg->next;
  119.     }
  120.     if ( !old.subpath_open && old.position_valid )
  121.         gx_path_add_point(ppath, old.position.x, old.position.y);
  122. #ifdef DEBUG
  123. if ( gs_debug['p'] )
  124.     gx_dump_path(ppath, "after copy_path");
  125. #endif
  126.     return 0;
  127. }
  128. /* Internal routine to flatten a curve. */
  129. /* This calls itself recursively, using binary subdivision, */
  130. /* until the approximation is good enough to satisfy the */
  131. /* flatness requirement.  The starting point is ppath->position, */
  132. /* which gets updated as line segments are added. */
  133.  
  134. /* Maximum number of points for sampling if we want accurate rasterizing. */
  135. /* (num_sample_max - 1)^3 must fit into an int. */
  136. #define num_sample_max (1 << ((sizeof(int) * 8 - 1) / 3))
  137.  
  138. /* Table of f(i) = 256 * sqrt(1 + (i/64)^2). */
  139. /* This is good to within 1% or better. */
  140. #define sqrt_index_shift 6        /* scaling of index */
  141. #define sqrt_value_shift 8        /* scaling of value */
  142. private int scaled_sqrt_tab[65] =
  143.    {    256, 256, 256, 256, 256, 256, 257, 257,
  144.     257, 258, 259, 259, 260, 261, 262, 262,
  145.     263, 264, 265, 267, 268, 269, 270, 272,
  146.     273, 274, 276, 277, 279, 281, 282, 284,
  147.     286, 288, 289, 291, 293, 295, 297, 299,
  148.     301, 304, 306, 308, 310, 312, 315, 317,
  149.     320, 322, 324, 327, 329, 332, 334, 337,
  150.     340, 342, 345, 348, 350, 353, 356, 359,
  151.     362
  152.    };
  153.  
  154. private int
  155. flatten_recur(gx_path *ppath,
  156.   fixed x1, fixed y1, fixed x2, fixed y2, fixed x3, fixed y3,
  157.   fixed scaled_flat)
  158. {    fixed
  159.       x0 = ppath->position.x,
  160.       y0 = ppath->position.y;
  161. top:
  162. #ifdef DEBUG
  163. if ( gs_debug['2'] )
  164.     dprintf4("[2]x0=%f y0=%f x1=%f y1=%f\n",
  165.          fixed2float(x0), fixed2float(y0),
  166.          fixed2float(x1), fixed2float(y1)),
  167.     dprintf4("   x2=%f y2=%f x3=%f y3=%f\n",
  168.          fixed2float(x2), fixed2float(y2),
  169.          fixed2float(x3), fixed2float(y3));
  170. #endif
  171.     /*
  172.      * Compute the maximum distance of the curve from
  173.      * the line (x0,y0)->(x3,y3).  We do this conservatively
  174.      * by observing that the curve is enclosed by the
  175.      * quadrilateral of its control points, so we simply
  176.      * compute the distances of (x1,y1) and (x2,y2)
  177.      * from the line.  Letting dx = x3-x0 and dy = y3-y0,
  178.      * the distance of (xp,yp) from the line is
  179.      * abs(N)/sqrt(D), where N = dy*(xp-x0)-dx*(yp-y0) and
  180.      * D = dx*dx+dy*dy; hence we want to test abs(N) <= sqrt(D)*F,
  181.      * where F is the flatness parameter from the graphics state.
  182.      * We can do this more efficiently by letting t=dy/dx, and
  183.      * testing abs(N1) <= sqrt(D1)*f, where N1=t*(xp-x0)-(yp-y0) and
  184.      * D1 = 1+t*t.  If dx < dy, we swap x and y for this
  185.      * computation.  This guarantees abs(t) <= 1, which allows us to
  186.      * compute sqrt(1+t*t) by table lookup on the high bits of abs(t).
  187.      *
  188.      * To avoid replacing shallow curves by short straight lines,
  189.      * we halve the flatness if the curve is very short.
  190.      * We don't do anything about long, nearly flat curves.
  191.      *
  192.      * Note that if scaled_flat is very small, we don't do any of this;
  193.      * instead, we just check whether abs(dx) and abs(dy) are
  194.      * small enough to switch over to sampling rather than dividing.
  195.       */
  196.      { fixed dx3 = x3 - x0;
  197.        fixed adx3 = any_abs(dx3);
  198.        fixed dy3 = y3 - y0;
  199.        fixed ady3 = any_abs(dy3);
  200.        fixed flat = scaled_flat;
  201.        /* We have to be quite careful to ensure that */
  202.        /* none of the multiplications will overflow. */
  203. #define short_max 0x7ff0L
  204. #define reduce_3(ad3, maxv)\
  205.   while ( ad3 > maxv )\
  206.     adx3 >>= 1, ady3 >>= 1,\
  207.     dx3 = arith_rshift_1(dx3), dy3 = arith_rshift_1(dy3)
  208. #define reduce_d(d)\
  209.   for ( shift = 0; (d < 0 ? d < -short_max : d > short_max); shift++ )\
  210.     d = arith_rshift_1(d)
  211.        if ( adx3 > ady3 )
  212.         {    fixed d, dx, dy, dist;
  213.         int shift;
  214.         if ( scaled_flat == 0 )
  215.         { if ( adx3 < int2fixed(num_sample_max / 2) )
  216.           { int n = fixed2int_var_rounded(adx3) * 2 + 1;
  217.             return flatten_sample(ppath, max(n, 3), x1, y1,
  218.                       x2, y2, x3, y3);
  219.           }
  220.           else goto sub;
  221.         }
  222.         if ( adx3 < 16 ) flat >>= 1;
  223.         reduce_3(ady3, short_max);
  224.         d = (scaled_sqrt_tab[(ady3 << sqrt_index_shift) / adx3] * flat) >> sqrt_value_shift;
  225.         dx = x1 - x0, dy = y1 - y0;
  226.         reduce_d(dx);
  227.         if ( ((dist = ((dx * dy3 / dx3) << shift) - dy) < 0 ?
  228.               -dist : dist) > d )
  229.           goto sub;    /* not flat enough */
  230.         dx = x2 - x0, dy = y2 - y0;
  231.         reduce_d(dx);
  232.         if ( ((dist = ((dx * dy3 / dx3) << shift) - dy) < 0 ?
  233.               -dist : dist) > d )
  234.           goto sub;    /* not flat enough */
  235.         }
  236.        else if ( ady3 != 0 )
  237.         {    fixed d, dy, dx, dist;
  238.         int shift;
  239.         if ( scaled_flat == 0 )
  240.         { if ( ady3 < int2fixed(num_sample_max / 2) )
  241.           { int n = fixed2int_var_rounded(ady3) * 2 + 1;
  242.             return flatten_sample(ppath, max(n, 3), x1, y1,
  243.                       x2, y2, x3, y3);
  244.           }
  245.           else goto sub;
  246.         }
  247.         if ( ady3 < 16 ) flat >>= 1;
  248.         reduce_3(adx3, short_max);
  249.         d = (scaled_sqrt_tab[(adx3 << sqrt_index_shift) / ady3] * flat) >> sqrt_value_shift;
  250.         dy = y1 - y0, dx = x1 - x0;
  251.         reduce_d(dy);
  252.         if ( ((dist = ((dy * dx3 / dy3) << shift) - dx) < 0 ?
  253.               -dist : dist) > d )
  254.           goto sub;    /* not flat enough */
  255.         dy = y2 - y0, dx = x2 - x0;
  256.         reduce_d(dy);
  257.         if ( ((dist = ((dy * dx3 / dy3) << shift) - dx) < 0 ?
  258.               -dist : dist) > d )
  259.           goto sub;    /* not flat enough */
  260.         }
  261.        else                /* adx3 = ady3 = 0 */
  262.         {    /* (x0,y0) is the same point as (x3,y3). */
  263.         /* This is an anomalous case.  If the entire curve */
  264.         /* is a single point, stop now, otherwise subdivide. */
  265.         if ( x1 != x0 || y1 != y0 || x2 != x0 || y2 != y0 )
  266.           goto sub;
  267.         }
  268.      }
  269.     /* Curve is flat enough.  Add a line and exit. */
  270.     if_debug2('2', "[2]\t*** x=%f, y=%f ***\n",
  271.           fixed2float(x3), fixed2float(y3));
  272.     return gx_path_add_line(ppath, x3, y3);
  273.  
  274.     /* Curve isn't flat enough.  Break into two pieces and recur. */
  275.     /* Algorithm is from "The Beta2-split: A special case of the */
  276.     /* Beta-spline Curve and Surface Representation," B. A. Barsky */
  277.     /* and A. D. DeRose, IEEE, 1985, courtesy of Crispin Goswell. */
  278. sub:
  279.     /* We have to define midpoint carefully to avoid overflow. */
  280.     /* (If it overflows, something really pathological is going on, */
  281.     /* but we could get infinite recursion that way....) */
  282. #define midpoint(a,b)\
  283.   (arith_rshift_1(a) + arith_rshift_1(b) + ((a) & (b) & 1))
  284.    {    fixed x01 = midpoint(x0, x1), y01 = midpoint(y0, y1);
  285.     fixed x12 = midpoint(x1, x2), y12 = midpoint(y1, y2);
  286.     fixed x02 = midpoint(x01, x12), y02 = midpoint(y01, y12);
  287.     int code;
  288.     /* Update x/y1, x/y2, and x/y0 now for the second half. */
  289.     x2 = midpoint(x2, x3), y2 = midpoint(y2, y3);
  290.     x1 = midpoint(x12, x2), y1 = midpoint(y12, y2);
  291.     code = flatten_recur(ppath, x01, y01, x02, y02,
  292.         (x0 = midpoint(x02, x1)), (y0 = midpoint(y02, y1)),
  293.         scaled_flat);
  294.     if ( code < 0 ) return code;
  295.    }    goto top;
  296. }
  297.  
  298. /* Flatten a segment of the path by repeated sampling. */
  299. /* n is the number of points to sample, including the endpoints; */
  300. /* we require n >= 3. */
  301. private int
  302. flatten_sample(gx_path *ppath, int n,
  303.   fixed x1, fixed y1, fixed x2, fixed y2, fixed x3, fixed y3)
  304. {    const fixed
  305.         x0 = ppath->position.x,
  306.         y0 = ppath->position.y;
  307.     /* We spell out some multiplies by 3, */
  308.     /* for the benefit of compilers that don't optimize this. */
  309.     const fixed
  310.         x01 = x1 - x0,
  311.         cx = (x01 << 1) + x01,        /* 3*(x1-x0) */
  312.         x12 = x2 - x1,
  313.         bx = (x12 << 1) + x12 - cx,    /* 3*(x2-2*x1+x0) */
  314.         ax = x3 - bx - cx - x0;        /* x3-3*x2+3*x1-x0 */
  315.     const fixed
  316.         y01 = y1 - y0,
  317.         cy = (y01 << 1) + y01,
  318.         y12 = y2 - y1,
  319.         by = (y12 << 1) + y12 - cy,
  320.         ay = y3 - by - cy - y0;
  321.     const int
  322.         n1 = n - 1,
  323.         n12 = n1 * n1,
  324.         n13 = n12 * n1;
  325.     fixed ptx = x0, pty = y0;
  326.     fixed x, y;
  327.     /*
  328.      * If all the coefficients lie between min_fast and max_fast,
  329.      * we can do everything in fixed point.  In this case we compute
  330.      * successive values by finite differences, using the formulas:
  331.         x(t) =
  332.           a*t^3 + b*t^2 + c*t + d =>
  333.         dx(t) = x(t+e)-x(t) =
  334.           a*(3*t^2*e + 3*t*e^2 + e^3) + b*(2*t*e + e^2) + c*e =
  335.           (3*a*e)*t^2 + (3*a*e^2 + 2*b*e)*t + (a*e^3 + b*e^2 + c*e) =>
  336.         d2x(t) = dx(t+e)-dx(t) =
  337.           (3*a*e)*(2*t*e + e^2) + (3*a*e^2 + 2*b*e)*e =
  338.           (6*a*e^2)*t + (6*a*e^3 + 2*b*e^2) =>
  339.         d3x(t) = d2x(t+e)-d2x(t) =
  340.           6*a*e^3;
  341.         x(0) = d, dx(0) = (a*e^3 + b*e^2 + c*e),
  342.           d2x(0) = 6*a*e^3 + 2*b*e^2;
  343.      * In these formulas, e = 1/n1; of course, there are separate
  344.      * computations for the x and y values.
  345.      */
  346. #define max_fast (max_fixed / 6)
  347. #define min_fast (-max_fast)
  348.     int fast;
  349.     float dt;        /* only if !fast */
  350.     int i;
  351.     /*
  352.      * We do exact rational arithmetic to avoid accumulating error.
  353.      * Each quantity is represented as I+R/n13, where I is an "integer"
  354.      * and the "remainder" R lies in the range 0 <= R < n13.  Note that
  355.      * R may temporarily exceed n13, and hence possibly overflow.
  356.      */
  357.     fixed idx, idy, id2x, id2y, id3x, id3y;        /* I */
  358.     int rx, ry, rdx, rdy, rd2x, rd2y, rd3x, rd3y;
  359.  
  360.     if_debug6('2', "[2]ax=%f bx=%f cx=%f\n   ay=%f by=%f cy=%f\n",
  361.           fixed2float(ax), fixed2float(bx), fixed2float(cx),
  362.           fixed2float(ay), fixed2float(by), fixed2float(cy));
  363. #define in_range(v) (v < max_fast && v > min_fast)
  364.     if ( n1 < num_sample_max &&    /* so n13 fits */
  365.          in_range(ax) && in_range(bx) && in_range(cx) &&
  366.          in_range(ay) && in_range(by) && in_range(cy)
  367.        )
  368.     {    fast = 1;
  369.         x = x0, y = y0;
  370.         rx = ry = 0;
  371.         /* Fast check for n == 3, a common special case */
  372.         /* for small characters. */
  373.         if ( n == 3 )
  374. #define poly2(a,b,c)\
  375.   arith_rshift_1(arith_rshift_1(arith_rshift_1(a) + b) + c)
  376.             idx = poly2(ax, bx, cx),
  377.             idy = poly2(ay, by, cy),
  378.             rdx = rdy = 0;
  379. #undef poly2
  380.         else
  381.         {    fixed bx2 = bx << 1, by2 = by << 1;
  382.             fixed ax6 = ((ax << 1) + ax) << 1,
  383.                   ay6 = ((ay << 1) + ay) << 1;
  384.             fixed qx, qy;
  385. #define adjust_rem(r, q)\
  386.   if ( r < 0 ) r += n13, q--
  387. #define adjust_rem_loop(r, q)\
  388.   while ( r < 0 ) r += n13, q--;\
  389.   while ( r >= n13 ) r -= n13, q++
  390.             /* We can compute all the remainders as ints, */
  391.             /* because we know they are less than n13. */
  392.             /* bx/y terms */
  393.             id2x = bx2 / n12, id2y = by2 / n12;
  394.             rd2x = ((int)bx2 - (int)id2x * n12) * n1,
  395.               rd2y = ((int)by2 - (int)id2y * n12) * n1;
  396.             idx = bx / n12, idy = by / n12;
  397.             rdx = ((int)bx - (int)idx * n12) * n1,
  398.               rdy = ((int)by - (int)idy * n12) * n1;
  399.             /* cx/y terms */
  400.             idx += qx = cx / n1, idy += qy = cy / n1;
  401.             rdx += ((int)cx - (int)qx * n1) * n12,
  402.               rdy += ((int)cy - (int)qy * n1) * n12;
  403.             /* ax/y terms */
  404.             idx += qx = ax / n13, idy += qy = ay / n13;
  405.             rdx += (int)ax - (int)qx * n13,
  406.               rdy += (int)ay - (int)qy * n13;
  407.             id2x += id3x = ax6 / n13, id2y += id3y = ay6 / n13;
  408.             rd2x += rd3x = (int)ax6 - (int)id3x * n13,
  409.               rd2y += rd3y = (int)ay6 - (int)id3y * n13;
  410.             adjust_rem_loop(rdx, idx);
  411.             adjust_rem_loop(rdy, idy);
  412.             adjust_rem_loop(rd2x, id2x);
  413.             adjust_rem_loop(rd2y, id2y);
  414.             adjust_rem(rd3x, id3x);
  415.             adjust_rem(rd3y, id3y);
  416. #undef adjust_rem
  417. #undef adjust_rem_loop
  418.         }
  419.     }
  420.     else
  421.         fast = 0, dt = 1.0 / (float)n1;
  422.     if_debug4('2', "[2]sampling %s n=%d\n[2]x=%g, y=%g\n",
  423.           (fast ? "fast" : "slow"), n,
  424.           fixed2float(x), fixed2float(y));
  425.     for ( i = 1; i < n1; i++ )
  426.     {    int code;
  427.         if ( fast )
  428.         {
  429. #ifdef DEBUG
  430. if ( gs_debug['2'] )
  431.             dprintf4("[2]dx=%f+%d, dy=%f+%d\n",
  432.                  fixed2float(idx), rdx,
  433.                  fixed2float(idy), rdy),
  434.             dprintf4("   d2x=%f+%d, d2y=%f+%d\n",
  435.                  fixed2float(id2x), rd2x,
  436.                  fixed2float(id2y), rd2y),
  437.             dprintf4("   d3x=%f+%d, d3y=%f+%d\n",
  438.                  fixed2float(id3x), rd3x,
  439.                  fixed2float(id3y), rd3y);
  440. #endif
  441. #define accum(i, r, di, dr)\
  442.   if ( (unsigned)(r += dr) >= (unsigned)n13 ) r -= n13, i++;\
  443.   i += di
  444.             accum(x, rx, idx, rdx);
  445.             accum(idx, rdx, id2x, rd2x);
  446.             accum(id2x, rd2x, id3x, rd3x);
  447.             accum(y, ry, idy, rdy);
  448.             accum(idy, rdy, id2y, rd2y);
  449.             accum(id2y, rd2y, id3y, rd3y);
  450. #undef accum
  451.         }
  452.         else
  453.         {    const float t = dt * (float)(i);
  454.             x = ((ax*t + bx)*t + cx)*t + x0;
  455.             y = ((ay*t + by)*t + cy)*t + y0;
  456.         }
  457.         if_debug3('2', "[2]%s x=%g, y=%g\n",
  458.               (((x ^ ptx) | (y ^ pty)) & float2fixed(-0.5) ?
  459.                "add" : "skip"),
  460.               fixed2float(x), fixed2float(y));
  461.         /* Skip very short segments */
  462.         if ( ((x ^ ptx) | (y ^ pty)) & float2fixed(-0.5) )
  463.         {    if ( (code = gx_path_add_line(ppath, x, y)) < 0 )
  464.                 return code;
  465.             ptx = x, pty = y;
  466.         }
  467.     }
  468.     if_debug2('2', "[2]last x=%g, y=%g\n",
  469.           fixed2float(x3), fixed2float(y3));
  470.     return gx_path_add_line(ppath, x3, y3);
  471. }
  472.  
  473. /*
  474.  *    The rest of this file is an analysis that will eventually
  475.  *    allow us to rasterize curves on the fly, by finding points
  476.  *    where Y reaches a local maximum or minimum, which allows us to
  477.  *    divide the curve into locally Y-monotonic sections.
  478.  */
  479.  
  480. /*
  481.     Let y(t) = a*t^3 + b*t^2 + c*t + d, 0 <= t <= 1.
  482.     Then dy(t) = 3*a*t^2 + 2*b*t + c.
  483.     y(t) has a local minimum or maximum (or inflection point)
  484.     precisely where dy(t) = 0.  Now the roots of dy(t) are
  485.         ( -2*b +/- sqrt(4*b^2 - 12*a*c) ) / 6*a
  486.        =    ( -b +/- sqrt(b*2 - 3*a*c) ) / 3*a
  487.     (Note that real roots exist iff b^2 >= 3*a*c.)
  488.     We want to know if these lie in the range (0..1).
  489.     (We don't care about the endpoints.)  Call such a root
  490.     a "valid zero."  We proceed as follows:
  491.         If sign(3*a + 2*b + c) ~= sign(c), a valid zero exists.
  492.         If sign(a) = sign(b), no valid zero exists.
  493.     Otherwise, we look for a local extremum of dy(t) by observing
  494.         d2y(t) = 6*a*t + 2*b
  495.     which has a zero only at
  496.         t1 = -b / 3*a
  497.     Now if t1 <= 0 or t1 >= 1, no valid zero exists.  Otherwise,
  498.     we compute
  499.         dy(t1) = c - (b^2 / 3*a)
  500.     Then a valid zero exists (at t1) iff sign(dy(t1)) ~= sign(c).
  501.  */
  502.