home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Graphics / Graphics.zip / povsrc31.zip / polysolv.c < prev    next >
C/C++ Source or Header  |  1999-10-20  |  32KB  |  1,712 lines
  1. /****************************************************************************
  2. *                polysolv.c
  3. *
  4. *  This file was written by Alexander Enzmann.  He wrote the code for
  5. *  4th-6th order shapes and generously provided us these enhancements.
  6. *
  7. *  from Persistence of Vision(tm) Ray Tracer
  8. *  Copyright 1996,1999 Persistence of Vision Team
  9. *---------------------------------------------------------------------------
  10. *  NOTICE: This source code file is provided so that users may experiment
  11. *  with enhancements to POV-Ray and to port the software to platforms other
  12. *  than those supported by the POV-Ray Team.  There are strict rules under
  13. *  which you are permitted to use this file.  The rules are in the file
  14. *  named POVLEGAL.DOC which should be distributed with this file.
  15. *  If POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  16. *  Team Coordinator by email to team-coord@povray.org or visit us on the web at
  17. *  http://www.povray.org. The latest version of POV-Ray may be found at this site.
  18. *
  19. * This program is based on the popular DKB raytracer version 2.12.
  20. * DKBTrace was originally written by David K. Buck.
  21. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  22. *
  23. *****************************************************************************/
  24.  
  25. #include "frame.h"
  26. #include "povray.h"
  27. #include "povproto.h"
  28. #include "vector.h"
  29. #include "polysolv.h"
  30.  
  31.  
  32.  
  33. /*****************************************************************************
  34. * Local preprocessor defines
  35. ******************************************************************************/
  36.  
  37. /*                  WARNING     WARNING    WARNING
  38.  
  39.    The following three constants have been defined so that quartic equations
  40.    will properly render on fpu/compiler combinations that only have 64 bit
  41.    IEEE precision.  Do not arbitrarily change any of these values.
  42.  
  43.    If you have a machine with a proper fpu/compiler combo - like a Mac with
  44.    a 68881, then use the native floating format (96 bits) and tune the
  45.    values for a little less tolerance: something like: factor1 = 1.0e15,
  46.    factor2 = -1.0e-7, factor3 = 1.0e-10.
  47.  
  48.    The meaning of the three constants are:
  49.       factor1 - the magnitude of difference between coefficients in a
  50.                 quartic equation at which the Sturmian root solver will
  51.                 kick in.  The Sturm solver is quite a bit slower than
  52.                 the algebraic solver, so the bigger this is, the faster
  53.                 the root solving will go.  The algebraic solver is less
  54.                 accurate so the bigger this is, the more likely you will
  55.                 get bad roots.
  56.  
  57.       factor2 - Tolerance value that defines how close the quartic equation
  58.                 is to being a square of a quadratic.  The closer this can
  59.                 get to zero before roots disappear, the less the chance
  60.                 you will get spurious roots.
  61.  
  62.       factor3 - Similar to factor2 at a later stage of the algebraic solver.
  63. */
  64.  
  65. #ifndef FUDGE_FACTOR1
  66.  #define FUDGE_FACTOR1 1.0e12
  67. #endif
  68. #ifndef FUDGE_FACTOR2
  69.  #define FUDGE_FACTOR2 -1.0e-5
  70. #endif
  71. #ifndef FUDGE_FACTOR3 
  72.  #define FUDGE_FACTOR3 1.0e-7
  73. #endif
  74.  
  75. /* Constants. */
  76.  
  77. #define TWO_M_PI_3  2.0943951023931954923084
  78. #define FOUR_M_PI_3 4.1887902047863909846168
  79.  
  80. /* Max number of iterations. */
  81.  
  82. #define MAX_ITERATIONS 50
  83.  
  84. /* A coefficient smaller than SMALL_ENOUGH is considered to be zero (0.0). */
  85.  
  86. #define SMALL_ENOUGH 1.0e-10
  87.  
  88. /* Smallest relative error we want. */
  89.  
  90. #define RELERROR 1.0e-12
  91.  
  92.  
  93.  
  94. /*****************************************************************************
  95. * Local typedefs
  96. ******************************************************************************/
  97.  
  98. typedef struct p
  99. {
  100.   int ord;
  101.   DBL coef[MAX_ORDER+1];
  102. }
  103. polynomial;
  104.  
  105.  
  106.  
  107. /*****************************************************************************
  108. * Local variables
  109. ******************************************************************************/
  110.  
  111.  
  112.  
  113. /*****************************************************************************
  114. * Static functions
  115. ******************************************************************************/
  116.  
  117. static int solve_quadratic (DBL *x, DBL *y);
  118. static int solve_cubic (DBL *x, DBL *y);
  119. static int solve_quartic (DBL *x, DBL *y);
  120. static int polysolve (int order, DBL *Coeffs, DBL *roots);
  121. static int modp (polynomial *u, polynomial *v, polynomial *r);
  122. static int regula_falsa (int order, DBL *coef, DBL a, DBL b, DBL *val);
  123. static int sbisect (int np, polynomial *sseq, DBL min, DBL max, int atmin, int atmax, DBL *roots);
  124. static int numchanges (int np, polynomial *sseq, DBL a);
  125. static DBL polyeval (DBL x, int n, DBL *Coeffs);
  126. static int buildsturm (int ord, polynomial *sseq);
  127. static int visible_roots (int np, polynomial *sseq, int *atneg, int *atpos);
  128. static int difficult_coeffs (int n, DBL *x);
  129.  
  130.  
  131.  
  132. /*****************************************************************************
  133. *
  134. * FUNCTION
  135. *
  136. *   modp
  137. *
  138. * INPUT
  139. *   
  140. * OUTPUT
  141. *   
  142. * RETURNS
  143. *   
  144. * AUTHOR
  145. *
  146. *   Alexander Enzmann
  147. *   
  148. * DESCRIPTION
  149. *
  150. *   Calculates the modulus of u(x) / v(x) leaving it in r.
  151. *   It returns 0 if r(x) is a constant.
  152. *
  153. *   NOTE: This function assumes the leading coefficient of v is 1 or -1.
  154. *
  155. * CHANGES
  156. *
  157. *   Okt 1996 : Added bug fix by Heiko Eissfeldt. [DB]
  158. *
  159. ******************************************************************************/
  160.  
  161. static int modp(polynomial *u, polynomial  *v, polynomial  *r)
  162. {
  163.   int k, j;
  164.  
  165.   *r = *u;
  166.  
  167.   if (v->coef[v->ord] < 0.0)
  168.   {
  169.     for (k = u->ord - v->ord - 1; k >= 0; k -= 2)
  170.     {
  171.       r->coef[k] = -r->coef[k];
  172.     }
  173.  
  174.     for (k = u->ord - v->ord; k >= 0; k--)
  175.     {
  176.       for (j = v->ord + k - 1; j >= k; j--)
  177.       {
  178.         r->coef[j] = -r->coef[j] - r->coef[v->ord + k] * v->coef[j - k];
  179.       }
  180.     }
  181.   }
  182.   else
  183.   {
  184.     for (k = u->ord - v->ord; k >= 0; k--)
  185.     {
  186.       for (j = v->ord + k - 1; j >= k; j--)
  187.       {
  188.         r->coef[j] -= r->coef[v->ord + k] * v->coef[j - k];
  189.       }
  190.     }
  191.   }
  192.  
  193.   k = v->ord - 1;
  194.  
  195.   while (k >= 0 && fabs(r->coef[k]) < SMALL_ENOUGH)
  196.   {
  197.     r->coef[k] = 0.0;
  198.  
  199.     k--;
  200.   }
  201.  
  202.   r->ord = (k < 0) ? 0 : k;
  203.  
  204.   return(r->ord);
  205. }
  206.  
  207.  
  208.  
  209. /*****************************************************************************
  210. *
  211. * FUNCTION
  212. *
  213. *   buildsturm
  214. *
  215. * INPUT
  216. *   
  217. * OUTPUT
  218. *   
  219. * RETURNS
  220. *   
  221. * AUTHOR
  222. *
  223. *   Alexander Enzmann
  224. *   
  225. * DESCRIPTION
  226. *
  227. *   Build the sturmian sequence for a polynomial.
  228. *
  229. * CHANGES
  230. *
  231. *   -
  232. *
  233. ******************************************************************************/
  234.  
  235. static int buildsturm(int ord, polynomial *sseq)
  236. {
  237.   int i;
  238.   DBL f, *fp, *fc;
  239.   polynomial *sp;
  240.  
  241.   sseq[0].ord = ord;
  242.   sseq[1].ord = ord - 1;
  243.  
  244.   /* calculate the derivative and normalize the leading coefficient. */
  245.  
  246.   f = fabs(sseq[0].coef[ord] * ord);
  247.  
  248.   fp = sseq[1].coef;
  249.   fc = sseq[0].coef + 1;
  250.  
  251.   for (i = 1; i <= ord; i++)
  252.   {
  253.     *fp++ = *fc++ * i / f;
  254.   }
  255.  
  256.   /* construct the rest of the Sturm sequence */
  257.  
  258.   for (sp = sseq + 2; modp(sp - 2, sp - 1, sp); sp++)
  259.   {
  260.     /* reverse the sign and normalize */
  261.  
  262.     f = -fabs(sp->coef[sp->ord]);
  263.  
  264.     for (fp = &sp->coef[sp->ord]; fp >= sp->coef; fp--)
  265.     {
  266.       *fp /= f;
  267.     }
  268.   }
  269.  
  270.   /* reverse the sign */
  271.  
  272.   sp->coef[0] = -sp->coef[0];
  273.  
  274.   return(sp - sseq);
  275. }
  276.  
  277.  
  278.  
  279. /*****************************************************************************
  280. *
  281. * FUNCTION
  282. *
  283. *   visible_roots
  284. *
  285. * INPUT
  286. *   
  287. * OUTPUT
  288. *   
  289. * RETURNS
  290. *
  291. * AUTHOR
  292. *
  293. *   Alexander Enzmann
  294. *   
  295. * DESCRIPTION
  296. *
  297. *   Find out how many visible intersections there are.
  298. *
  299. * CHANGES
  300. *
  301. *   -
  302. *
  303. ******************************************************************************/
  304.  
  305. static int visible_roots(int np, polynomial *sseq, int *atzer, int  *atpos)
  306. {
  307.   int atposinf, atzero;
  308.   polynomial *s;
  309.   DBL f, lf;
  310.  
  311.   atposinf = atzero = 0;
  312.  
  313.   /* changes at positve infinity */
  314.  
  315.   lf = sseq[0].coef[sseq[0].ord];
  316.  
  317.   for (s = sseq + 1; s <= sseq + np; s++)
  318.   {
  319.     f = s->coef[s->ord];
  320.  
  321.     if (lf == 0.0 || lf * f < 0)
  322.     {
  323.       atposinf++;
  324.     }
  325.  
  326.     lf = f;
  327.   }
  328.  
  329.   /* Changes at zero */
  330.  
  331.   lf = sseq[0].coef[0];
  332.  
  333.   for (s = sseq+1; s <= sseq + np; s++)
  334.   {
  335.     f = s->coef[0];
  336.  
  337.     if (lf == 0.0 || lf * f < 0)
  338.     {
  339.       atzero++;
  340.     }
  341.  
  342.     lf = f;
  343.   }
  344.  
  345.   *atzer = atzero;
  346.   *atpos = atposinf;
  347.  
  348.   return(atzero - atposinf);
  349. }
  350.  
  351.  
  352.  
  353. /*****************************************************************************
  354. *
  355. * FUNCTION
  356. *
  357. *   numchanges
  358. *
  359. * INPUT
  360. *   
  361. * OUTPUT
  362. *   
  363. * RETURNS
  364. *   
  365. * AUTHOR
  366. *
  367. *   Alexander Enzmann
  368. *   
  369. * DESCRIPTION
  370. *
  371. *   Return the number of sign changes in the Sturm sequence in
  372. *   sseq at the value a.
  373. *
  374. * CHANGES
  375. *
  376. *   -
  377. *
  378. ******************************************************************************/
  379.  
  380. static int numchanges(int np, polynomial *sseq, DBL a)
  381. {
  382.   int changes;
  383.   DBL f, lf;
  384.   polynomial   *s;
  385.   changes = 0;
  386.  
  387.   lf = polyeval(a, sseq[0].ord, sseq[0].coef);
  388.  
  389.   for (s = sseq + 1; s <= sseq + np; s++)
  390.   {
  391.     f = polyeval(a, s->ord, s->coef);
  392.  
  393.     if (lf == 0.0 || lf * f < 0)
  394.     {
  395.       changes++;
  396.     }
  397.  
  398.     lf = f;
  399.   }
  400.  
  401.   return(changes);
  402. }
  403.  
  404.  
  405.  
  406. /*****************************************************************************
  407. *
  408. * FUNCTION
  409. *
  410. *   sbisect
  411. *
  412. * INPUT
  413. *   
  414. * OUTPUT
  415. *   
  416. * RETURNS
  417. *   
  418. * AUTHOR
  419. *
  420. *   Alexander Enzmann
  421. *   
  422. * DESCRIPTION
  423. *
  424. *   Uses a bisection based on the sturm sequence for the polynomial
  425. *   described in sseq to isolate intervals in which roots occur,
  426. *   the roots are returned in the roots array in order of magnitude.
  427. *
  428. *   NOTE: This routine has one severe bug: When the interval containing the
  429. *         root [min, max] has a root at one of its endpoints, as well as one
  430. *         within the interval, the root at the endpoint will be returned
  431. *         rather than the one inside.
  432. *
  433. * CHANGES
  434. *
  435. *   -
  436. *
  437. ******************************************************************************/
  438.  
  439. static int sbisect(int np, polynomial *sseq, DBL min_value, DBL  max_value, int atmin, int  atmax, DBL *roots)
  440. {
  441.   DBL mid;
  442.   int n1, n2, its, atmid;
  443.   
  444.   if ((atmin - atmax) == 1)
  445.   {
  446.     /* first try using regula-falsa to find the root.  */
  447.  
  448.     if (regula_falsa(sseq->ord, sseq->coef, min_value, max_value, roots))
  449.     {
  450.       return(1);
  451.     }
  452.     else
  453.     {
  454.       /* That failed, so now find it by bisection */
  455.  
  456.       for (its = 0; its < MAX_ITERATIONS; its++)
  457.       {
  458.         mid = (min_value + max_value) / 2;
  459.  
  460.         atmid = numchanges(np, sseq, mid);
  461.  
  462.         /* The follow only happens if there is a bug.  And
  463.            unfortunately, there is. CEY 04/97 
  464.          */
  465.         if ((atmid<atmax) || (atmid>atmin))
  466.         {
  467.            return(0);
  468.         }
  469.   
  470.         if (fabs(mid) > RELERROR)
  471.         {
  472.           if (fabs((max_value - min_value) / mid) < RELERROR)
  473.           {
  474.             roots[0] = mid;
  475.  
  476.             return(1);
  477.           }
  478.         }
  479.         else
  480.         {
  481.           if (fabs(max_value - min_value) < RELERROR)
  482.           {
  483.             roots[0] = mid;
  484.  
  485.             return(1);
  486.           }
  487.         }
  488.  
  489.         if ((atmin - atmid) == 0)
  490.         {
  491.           min_value = mid;
  492.         }
  493.         else
  494.         {
  495.           max_value = mid;
  496.         }
  497.       }
  498.  
  499.       /* Bisection took too long - just return what we got */
  500.  
  501.       roots[0] = mid;
  502.  
  503.       return(1);
  504.     }
  505.   }
  506.  
  507.   /* There is more than one root in the interval.
  508.      Bisect to find new intervals. */
  509.  
  510.   for (its = 0; its < MAX_ITERATIONS; its++)
  511.   {
  512.     mid = (min_value + max_value) / 2;
  513.  
  514.     atmid = numchanges(np, sseq, mid);
  515.  
  516.     /* The follow only happens if there is a bug.  And
  517.        unfortunately, there is. CEY 04/97 
  518.      */
  519.     if ((atmid<atmax) || (atmid>atmin))
  520.     {
  521.        return(0);
  522.     }
  523.  
  524.     if (fabs(mid) > RELERROR)
  525.     {
  526.       if (fabs((max_value - min_value) / mid) < RELERROR)
  527.       {
  528.         roots[0] = mid;
  529.  
  530.         return(1);
  531.       }
  532.     }
  533.     else
  534.     {
  535.       if (fabs(max_value - min_value) < RELERROR)
  536.       {
  537.         roots[0] = mid;
  538.  
  539.         return(1);
  540.       }
  541.     }
  542.   
  543.     n1 = atmin - atmid;
  544.     n2 = atmid - atmax;
  545.  
  546.     if ((n1 != 0) && (n2 != 0))
  547.     {
  548.       n1 = sbisect(np, sseq, min_value, mid, atmin, atmid, roots);
  549.       n2 = sbisect(np, sseq, mid, max_value, atmid, atmax, &roots[n1]);
  550.  
  551.       return(n1 + n2);
  552.     }
  553.  
  554.     if (n1 == 0)
  555.     {
  556.       min_value = mid;
  557.     }
  558.     else
  559.     {
  560.       max_value = mid;
  561.     }
  562.   }
  563.  
  564.   /* Took too long to bisect.  Just return what we got. */
  565.  
  566.   roots[0] = mid;
  567.  
  568.   return(1);
  569. }
  570.  
  571.  
  572.  
  573. /*****************************************************************************
  574. *
  575. * FUNCTION
  576. *
  577. *   polyeval
  578. *
  579. * INPUT
  580. *   
  581. * OUTPUT
  582. *   
  583. * RETURNS
  584. *   
  585. * AUTHOR
  586. *
  587. *   Alexander Enzmann
  588. *   
  589. * DESCRIPTION
  590. *
  591. *   Evaluate the value of a polynomial at the given value x.
  592. *
  593. *   The coefficients are stored in c in the following order:
  594. *
  595. *     c[0] + c[1] * x + c[2] * x ^ 2 + c[3] * x ^ 3 + ...
  596. *
  597. * CHANGES
  598. *
  599. *   -
  600. *
  601. ******************************************************************************/
  602.  
  603. static DBL polyeval(DBL x, int n, DBL  *Coeffs)
  604. {
  605.   register int i;
  606.   DBL val;
  607.  
  608.   val = Coeffs[n];
  609.  
  610.   for (i = n-1; i >= 0; i--)
  611.   {
  612.     val = val * x + Coeffs[i];
  613.   }
  614.  
  615.   return(val);
  616. }
  617.  
  618.  
  619.  
  620. /*****************************************************************************
  621. *
  622. * FUNCTION
  623. *
  624. *   regular_falsa
  625. *
  626. * INPUT
  627. *   
  628. * OUTPUT
  629. *
  630. * RETURNS
  631. *
  632. * AUTHOR
  633. *
  634. *   Alexander Enzmann
  635. *
  636. * DESCRIPTION
  637. *
  638. *   Close in on a root by using regula-falsa.
  639. *
  640. * CHANGES
  641. *
  642. *   -
  643. *
  644. ******************************************************************************/
  645.  
  646. static int regula_falsa(int order, DBL *coef, DBL a, DBL  b, DBL  *val)
  647. {
  648.   int its;
  649.   DBL fa, fb, x, fx, lfx;
  650.  
  651.   fa = polyeval(a, order, coef);
  652.   fb = polyeval(b, order, coef);
  653.  
  654.   if (fa * fb > 0.0)
  655.   {
  656.     return 0;
  657.   }
  658.  
  659.   if (fabs(fa) < SMALL_ENOUGH)
  660.   {
  661.     *val = a;
  662.  
  663.     return(1);
  664.   }
  665.  
  666.   if (fabs(fb) < SMALL_ENOUGH)
  667.   {
  668.     *val = b;
  669.  
  670.     return(1);
  671.   }
  672.  
  673.   lfx = fa;
  674.  
  675.   for (its = 0; its < MAX_ITERATIONS; its++)
  676.   {
  677.     x = (fb * a - fa * b) / (fb - fa);
  678.  
  679.     fx = polyeval(x, order, coef);
  680.  
  681.     if (fabs(x) > RELERROR)
  682.     {
  683.       if (fabs(fx / x) < RELERROR)
  684.       {
  685.         *val = x;
  686.  
  687.         return(1);
  688.       }
  689.     }
  690.     else
  691.     {
  692.       if (fabs(fx) < RELERROR)
  693.       {
  694.         *val = x;
  695.  
  696.         return(1);
  697.       }
  698.     }
  699.  
  700.     if (fa < 0)
  701.     {
  702.       if (fx < 0)
  703.       {
  704.         a = x;
  705.  
  706.         fa = fx;
  707.  
  708.         if ((lfx * fx) > 0)
  709.         {
  710.           fb /= 2;
  711.         }
  712.       }
  713.       else
  714.       {
  715.         b = x;
  716.  
  717.         fb = fx;
  718.  
  719.         if ((lfx * fx) > 0)
  720.         {
  721.           fa /= 2;
  722.         }
  723.       }
  724.     }
  725.     else
  726.     {
  727.       if (fx < 0)
  728.       {
  729.         b = x;
  730.  
  731.         fb = fx;
  732.  
  733.         if ((lfx * fx) > 0)
  734.         {
  735.           fa /= 2;
  736.         }
  737.       }
  738.       else
  739.       {
  740.         a = x;
  741.  
  742.         fa = fx;
  743.  
  744.         if ((lfx * fx) > 0)
  745.         {
  746.           fb /= 2;
  747.         }
  748.       }
  749.     }
  750.  
  751.     /* Check for underflow in the domain */
  752.  
  753.     if (fabs(b-a) < RELERROR)
  754.     {
  755.       *val = x;
  756.  
  757.       return(1);
  758.     }
  759.  
  760.     lfx = fx;
  761.   }
  762.  
  763.   return(0);
  764. }
  765.  
  766.  
  767.  
  768. /*****************************************************************************
  769. *
  770. * FUNCTION
  771. *
  772. *   solve_quadratic
  773. *
  774. * INPUT
  775. *
  776. * OUTPUT
  777. *
  778. * RETURNS
  779. *   
  780. * AUTHOR
  781. *
  782. *   Alexander Enzmann
  783. *   
  784. * DESCRIPTION
  785. *
  786. *   Solve the quadratic equation:
  787. *
  788. *     x[0] * x^2 + x[1] * x + x[2] = 0.
  789. *
  790. *   The value returned by this function is the number of real roots.
  791. *   The roots themselves are returned in y[0], y[1].
  792. *
  793. * CHANGES
  794. *
  795. *   -
  796. *
  797. ******************************************************************************/
  798.  
  799. static int solve_quadratic(DBL *x, DBL  *y)
  800. {
  801.   DBL d, t, a, b, c;
  802.  
  803.   a = x[0];
  804.   b = -x[1];
  805.   c = x[2];
  806.  
  807.   if (a == 0.0)
  808.   {
  809.     if (b == 0.0)
  810.     {
  811.       return(0);
  812.     }
  813.  
  814.     y[0] = c / b;
  815.  
  816.     return(1);
  817.   }
  818.  
  819.   d = b * b - 4.0 * a * c;
  820.  
  821.   /* Treat values of d around 0 as 0. */
  822.  
  823.   if ((d > -SMALL_ENOUGH) && (d < SMALL_ENOUGH))
  824.   {
  825.     y[0] = 0.5 * b / a;
  826.  
  827.     return(1);
  828.   }
  829.   else
  830.   {
  831.     if (d < 0.0)
  832.     {
  833.       return(0);
  834.     }
  835.   }
  836.  
  837.   d = sqrt(d);
  838.  
  839.   t = 2.0 * a;
  840.  
  841.   y[0] = (b + d) / t;
  842.   y[1] = (b - d) / t;
  843.  
  844.   return(2);
  845. }
  846.  
  847.  
  848.  
  849. /*****************************************************************************
  850. *
  851. * FUNCTION
  852. *
  853. *   solve_cubic
  854. *
  855. * INPUT
  856. *   
  857. * OUTPUT
  858. *   
  859. * RETURNS
  860. *   
  861. * AUTHOR
  862. *
  863. *   Alexander Enzmann
  864. *   
  865. * DESCRIPTION
  866. *
  867. *
  868. *   Solve the cubic equation:
  869. *
  870. *     x[0] * x^3 + x[1] * x^2 + x[2] * x + x[3] = 0.
  871. *
  872. *   The result of this function is an integer that tells how many real
  873. *   roots exist.  Determination of how many are distinct is up to the
  874. *   process that calls this routine.  The roots that exist are stored
  875. *   in (y[0], y[1], y[2]).
  876. *
  877. *   NOTE: This function relies very heavily on trigonometric functions and
  878. *         the square root function.  If an alternative solution is found
  879. *         that does not rely on transcendentals this code will be replaced.
  880. *
  881. * CHANGES
  882. *
  883. *   -
  884. *
  885. ******************************************************************************/
  886.  
  887. static int solve_cubic(DBL *x, DBL  *y)
  888. {
  889.   DBL Q, R, Q3, R2, sQ, d, an, theta;
  890.   DBL A2, a0, a1, a2, a3;
  891.  
  892.   a0 = x[0];
  893.  
  894.   if (a0 == 0.0)
  895.   {
  896.     return(solve_quadratic(&x[1], y));
  897.   }
  898.   else
  899.   {
  900.     if (a0 != 1.0)
  901.     {
  902.       a1 = x[1] / a0;
  903.       a2 = x[2] / a0;
  904.       a3 = x[3] / a0;
  905.     }
  906.     else
  907.     {
  908.       a1 = x[1];
  909.       a2 = x[2];
  910.       a3 = x[3];
  911.     }
  912.   }
  913.  
  914.   A2 = a1 * a1;
  915.  
  916.   Q = (A2 - 3.0 * a2) / 9.0;
  917.  
  918.   /* Modified to save some multiplications and to avoid a floating point
  919.      exception that occured with DJGPP and full optimization. [DB 8/94] */
  920.  
  921.   R = (a1 * (A2 - 4.5 * a2) + 13.5 * a3) / 27.0;
  922.  
  923.   Q3 = Q * Q * Q;
  924.  
  925.   R2 = R * R;
  926.  
  927.   d = Q3 - R2;
  928.  
  929.   an = a1 / 3.0;
  930.  
  931.   if (d >= 0.0)
  932.   {
  933.     /* Three real roots. */
  934.  
  935.     d = R / sqrt(Q3);
  936.  
  937.     theta = acos(d) / 3.0;
  938.  
  939.     sQ = -2.0 * sqrt(Q);
  940.  
  941.     y[0] = sQ * cos(theta) - an;
  942.     y[1] = sQ * cos(theta + TWO_M_PI_3) - an;
  943.     y[2] = sQ * cos(theta + FOUR_M_PI_3) - an;
  944.  
  945.     return(3);
  946.   }
  947.   else
  948.   {
  949.     sQ = pow(sqrt(R2 - Q3) + fabs(R), 1.0 / 3.0);
  950.  
  951.     if (R < 0)
  952.     {
  953.       y[0] = (sQ + Q / sQ) - an;
  954.     }
  955.     else
  956.     {
  957.       y[0] = -(sQ + Q / sQ) - an;
  958.     }
  959.  
  960.     return(1);
  961.   }
  962. }
  963.  
  964. #ifdef USE_NEW_DIFFICULT_COEFFS
  965.  
  966. /*****************************************************************************
  967. *
  968. * FUNCTION
  969. *
  970. *   difficult_coeffs
  971. *
  972. * INPUT
  973. *   
  974. * OUTPUT
  975. *   
  976. * RETURNS
  977. *   
  978. * AUTHOR
  979. *
  980. *   Alexander Enzmann
  981. *
  982. * DESCRIPTION
  983. *
  984. *   Test to see if any coeffs are more than 6 orders of magnitude
  985. *   larger than the smallest.
  986. *
  987. * CHANGES
  988. *
  989. *   -
  990. *
  991. ******************************************************************************/
  992.  
  993. static int difficult_coeffs(int n, DBL *x)
  994. {
  995.   int i, flag = 0;
  996.   DBL t, biggest;
  997.  
  998.   biggest = fabs(x[0]);
  999.  
  1000.   for (i = 1; i <= n; i++)
  1001.   {
  1002.     t = fabs(x[i]);
  1003.  
  1004.     if (t > biggest)
  1005.     {
  1006.       biggest = t;
  1007.     }
  1008.   }
  1009.  
  1010.   /* Everything is zero no sense in doing any more */
  1011.  
  1012.   if (biggest == 0.0)
  1013.   {
  1014.     return(0);
  1015.   }
  1016.  
  1017.   for (i = 0; i <= n; i++)
  1018.   {
  1019.     if (x[i] != 0.0)
  1020.     {
  1021.       if (fabs(biggest / x[i]) > FUDGE_FACTOR1)
  1022.       {
  1023.         x[i] = 0.0;
  1024.         flag = 1;
  1025.       }
  1026.     }
  1027.   }
  1028.  
  1029.   return(flag);
  1030. }
  1031. #else
  1032. /*****************************************************************************
  1033. *
  1034. * FUNCTION
  1035. *
  1036. *   difficult_coeffs
  1037. *
  1038. * INPUT
  1039. *   
  1040. * OUTPUT
  1041. *   
  1042. * RETURNS
  1043. *   
  1044. * AUTHOR
  1045. *
  1046. *   Alexander Enzmann
  1047. *   
  1048. * DESCRIPTION
  1049. *
  1050. *   Test to see if any coeffs are more than 6 orders of magnitude
  1051. *   larger than the smallest (function from POV 1.0) [DB 8/94].
  1052. *
  1053. * CHANGES
  1054. *
  1055. *   -
  1056. *
  1057. ******************************************************************************/
  1058.  
  1059. static int difficult_coeffs(int n, DBL *x)
  1060. {
  1061.   int i;
  1062.   DBL biggest;
  1063.  
  1064.   biggest = 0.0;
  1065.  
  1066.   for (i = 0; i <= n; i++)
  1067.   {
  1068.     if (fabs(x[i]) > biggest)
  1069.     {
  1070.       biggest = x[i];
  1071.     }
  1072.   }
  1073.  
  1074.   /* Everything is zero no sense in doing any more */
  1075.  
  1076.   if (biggest == 0.0)
  1077.   {
  1078.     return(0);
  1079.   }
  1080.  
  1081.   for (i = 0; i <= n; i++)
  1082.   {
  1083.     if (x[i] != 0.0)
  1084.     {
  1085.       if (fabs(biggest / x[i]) > FUDGE_FACTOR1)
  1086.       {
  1087.         return(1);
  1088.       }
  1089.     }
  1090.   }
  1091.  
  1092.   return(0);
  1093. }
  1094.  
  1095. #endif
  1096.  
  1097. #ifdef TEST_SOLVER
  1098. /*****************************************************************************
  1099. *
  1100. * FUNCTION
  1101. *
  1102. *   solve_quartic
  1103. *
  1104. * INPUT
  1105. *   
  1106. * OUTPUT
  1107. *   
  1108. * RETURNS
  1109. *   
  1110. * AUTHOR
  1111. *
  1112. *   Alexander Enzmann
  1113. *
  1114. * DESCRIPTION
  1115. *
  1116. *   The old way of solving quartics algebraically.
  1117. *   This is an adaptation of the method of Lodovico Ferrari (Circa 1731).
  1118. *
  1119. * CHANGES
  1120. *
  1121. *   -
  1122. *
  1123. ******************************************************************************/
  1124.  
  1125. static int solve_quartic(DBL *x, DBL  *results)
  1126. {
  1127.   DBL cubic[4], roots[3];
  1128.   DBL a0, a1, y, d1, x1, t1, t2;
  1129.   DBL c0, c1, c2, c3, c4, d2, q1, q2;
  1130.   int i;
  1131.  
  1132.   c0 = x[0];
  1133.  
  1134.   if (c0 != 1.0)
  1135.   {
  1136.     c1 = x[1] / c0;
  1137.     c2 = x[2] / c0;
  1138.     c3 = x[3] / c0;
  1139.     c4 = x[4] / c0;
  1140.   }
  1141.   else
  1142.   {
  1143.     c1 = x[1];
  1144.     c2 = x[2];
  1145.     c3 = x[3];
  1146.     c4 = x[4];
  1147.   }
  1148.  
  1149.   /* The first step is to take the original equation:
  1150.  
  1151.        x^4 + b*x^3 + c*x^2 + d*x + e = 0
  1152.  
  1153.      and rewrite it as:
  1154.  
  1155.        x^4 + b*x^3 = -c*x^2 - d*x - e,
  1156.  
  1157.      adding (b*x/2)^2 + (x^2 + b*x/2)y + y^2/4 to each side gives a
  1158.      perfect square on the lhs:
  1159.  
  1160.        (x^2 + b*x/2 + y/2)^2 = (b^2/4 - c + y)x^2 + (b*y/2 - d)x + y^2/4 - e
  1161.  
  1162.      By choosing the appropriate value for y, the rhs can be made a perfect
  1163.      square also.  This value is found when the rhs is treated as a quadratic
  1164.      in x with the discriminant equal to 0.  This will be true when:
  1165.  
  1166.        (b*y/2 - d)^2 - 4.0 * (b^2/4 - c*y)*(y^2/4 - e) = 0, or
  1167.  
  1168.        y^3 - c*y^2 + (b*d - 4*e)*y - b^2*e + 4*c*e - d^2 = 0.
  1169.  
  1170.      This is called the resolvent of the quartic equation.  */
  1171.  
  1172.   a0 = 4.0 * c4;
  1173.  
  1174.   cubic[0] = 1.0;
  1175.   cubic[1] = -1.0 * c2;
  1176.   cubic[2] = c1 * c3 - a0;
  1177.   cubic[3] = a0 * c2 - c1 * c1 * c4 - c3 * c3;
  1178.  
  1179.   i = solve_cubic(&cubic[0], &roots[0]);
  1180.  
  1181.   if (i > 0)
  1182.   {
  1183.     y = roots[0];
  1184.   }
  1185.   else
  1186.   {
  1187.     return(0);
  1188.   }
  1189.  
  1190.   /* What we are left with is a quadratic squared on the lhs and a
  1191.      linear term on the right.  The linear term has one of two signs,
  1192.      take each and add it to the lhs.  The form of the quartic is now:
  1193.  
  1194.        a' = b^2/4 - c + y,    b' = b*y/2 - d, (from rhs quadritic above)
  1195.  
  1196.        (x^2 + b*x/2 + y/2) = +sqrt(a'*(x + 1/2 * b'/a')^2), and
  1197.        (x^2 + b*x/2 + y/2) = -sqrt(a'*(x + 1/2 * b'/a')^2).
  1198.  
  1199.      By taking the linear term from each of the right hand sides and
  1200.      adding to the appropriate part of the left hand side, two quadratic
  1201.      formulas are created.  By solving each of these the four roots of
  1202.      the quartic are determined.
  1203.   */
  1204.  
  1205.   i = 0;
  1206.  
  1207.   a0 = c1 / 2.0;
  1208.   a1 = y / 2.0;
  1209.  
  1210.   t1 = a0 * a0 - c2 + y;
  1211.  
  1212.   if (t1 < 0.0)
  1213.   {
  1214.     if (t1 > FUDGE_FACTOR2)
  1215.     {
  1216.       t1 = 0.0;
  1217.     }
  1218.     else
  1219.     {
  1220.       /* First Special case, a' < 0 means all roots are complex. */
  1221.  
  1222.       return(0);
  1223.       }
  1224.     }
  1225.   }
  1226.  
  1227.   if (t1 < FUDGE_FACTOR3)
  1228. {
  1229.     /* Second special case, the "x" term on the right hand side above
  1230.        has vanished.  In this case:
  1231.  
  1232.          (x^2 + b*x/2 + y/2) = +sqrt(y^2/4 - e), and
  1233.          (x^2 + b*x/2 + y/2) = -sqrt(y^2/4 - e).  */
  1234.  
  1235.     t2 = a1 * a1 - c4;
  1236.  
  1237.     if (t2 < 0.0)
  1238.     {
  1239.       return(0);
  1240.     }
  1241.  
  1242.     x1 = 0.0;
  1243.     d1 = sqrt(t2);
  1244.   }
  1245.   else
  1246.   {
  1247.     x1 = sqrt(t1);
  1248.     d1 = 0.5 * (a0 * y - c3) / x1;
  1249.   }
  1250.  
  1251.   /* Solve the first quadratic */
  1252.  
  1253.   q1 = -a0 - x1;
  1254.   q2 = a1 + d1;
  1255.   d2 = q1 * q1 - 4.0 * q2;
  1256.  
  1257.   if (d2 >= 0.0)
  1258. {
  1259.     d2 = sqrt(d2);
  1260.  
  1261.     results[0] = 0.5 * (q1 + d2);
  1262.     results[1] = 0.5 * (q1 - d2);
  1263.  
  1264.     i = 2;
  1265.   }
  1266.  
  1267.   /* Solve the second quadratic */
  1268.  
  1269.   q1 = q1 + x1 + x1;
  1270.   q2 = a1 - d1;
  1271.   d2 = q1 * q1 - 4.0 * q2;
  1272.  
  1273.   if (d2 >= 0.0)
  1274. {
  1275.     d2 = sqrt(d2);
  1276.  
  1277.     results[i++] = 0.5 * (q1 + d2);
  1278.     results[i++] = 0.5 * (q1 - d2);
  1279.   }
  1280.  
  1281.   return(i);
  1282. }
  1283. #else
  1284. /*****************************************************************************
  1285. *
  1286. * FUNCTION
  1287. *
  1288. *   solve_quartic
  1289. *
  1290. * INPUT
  1291. *   
  1292. * OUTPUT
  1293. *   
  1294. * RETURNS
  1295. *   
  1296. * AUTHOR
  1297. *
  1298. *   Alexander Enzmann
  1299. *
  1300. * DESCRIPTION
  1301. *
  1302. *   Solve a quartic using the method of Francois Vieta (Circa 1735).
  1303. *
  1304. * CHANGES
  1305. *
  1306. *   -
  1307. *
  1308. ******************************************************************************/
  1309.  
  1310. static int solve_quartic(DBL *x, DBL *results)
  1311. {
  1312.   DBL cubic[4], roots[3];
  1313.   DBL c12, z, p, q, q1, q2, r, d1, d2;
  1314.   DBL c0, c1, c2, c3, c4;
  1315.   int i;
  1316.  
  1317.   /* Make sure the quartic has a leading coefficient of 1.0 */
  1318.  
  1319.   c0 = x[0];
  1320.  
  1321.   if (c0 != 1.0)
  1322. {
  1323.     c1 = x[1] / c0;
  1324.     c2 = x[2] / c0;
  1325.     c3 = x[3] / c0;
  1326.     c4 = x[4] / c0;
  1327.   }
  1328.   else
  1329.   {
  1330.     c1 = x[1];
  1331.     c2 = x[2];
  1332.     c3 = x[3];
  1333.     c4 = x[4];
  1334.   }
  1335.  
  1336.   /* Compute the cubic resolvant */
  1337.  
  1338.   c12 = c1 * c1;
  1339.   p = -0.375 * c12 + c2;
  1340.   q = 0.125 * c12 * c1 - 0.5 * c1 * c2 + c3;
  1341.   r = -0.01171875 * c12 * c12 + 0.0625 * c12 * c2 - 0.25 * c1 * c3 + c4;
  1342.  
  1343.   cubic[0] = 1.0;
  1344.   cubic[1] = -0.5 * p;
  1345.   cubic[2] = -r;
  1346.   cubic[3] = 0.5 * r * p - 0.125 * q * q;
  1347.  
  1348.   i = solve_cubic(cubic, roots);
  1349.  
  1350.   if (i > 0)
  1351. {
  1352.     z = roots[0];
  1353.   }
  1354.   else
  1355.   {
  1356.     return(0);
  1357.   }
  1358.  
  1359.   d1 = 2.0 * z - p;
  1360.  
  1361.   if (d1 < 0.0)
  1362. {
  1363.     if (d1 > -SMALL_ENOUGH)
  1364.     {
  1365.       d1 = 0.0;
  1366.     }
  1367.     else
  1368.     {
  1369.       return(0);
  1370.     }
  1371.   }
  1372.  
  1373.   if (d1 < SMALL_ENOUGH)
  1374. {
  1375.     d2 = z * z - r;
  1376.  
  1377.     if (d2 < 0.0)
  1378.     {
  1379.       return(0);
  1380.     }
  1381.  
  1382.     d2 = sqrt(d2);
  1383.   }
  1384.   else
  1385.   {
  1386.     d1 = sqrt(d1);
  1387.     d2 = 0.5 * q / d1;
  1388.   }
  1389.  
  1390.   /* Set up useful values for the quadratic factors */
  1391.  
  1392.   q1 = d1 * d1;
  1393.   q2 = -0.25 * c1;
  1394.  
  1395.   i = 0;
  1396.  
  1397.   /* Solve the first quadratic */
  1398.  
  1399.   p = q1 - 4.0 * (z - d2);
  1400.  
  1401.   if (p == 0)
  1402. {
  1403.     results[i++] = -0.5 * d1 - q2;
  1404.   }
  1405.   else
  1406.   {
  1407.     if (p > 0)
  1408.     {
  1409.       p = sqrt(p);
  1410.       results[i++] = -0.5 * (d1 + p) + q2;
  1411.       results[i++] = -0.5 * (d1 - p) + q2;
  1412.     }
  1413.   }
  1414.  
  1415.   /* Solve the second quadratic */
  1416.  
  1417.   p = q1 - 4.0 * (z + d2);
  1418.  
  1419.   if (p == 0)
  1420. {
  1421.     results[i++] = 0.5 * d1 - q2;
  1422.   }
  1423.   else
  1424.   {
  1425.     if (p > 0)
  1426.     {
  1427.       p = sqrt(p);
  1428.       results[i++] = 0.5 * (d1 + p) + q2;
  1429.       results[i++] = 0.5 * (d1 - p) + q2;
  1430.     }
  1431.   }
  1432.  
  1433.   return(i);
  1434. }
  1435. #endif
  1436.  
  1437.  
  1438.  
  1439. /*****************************************************************************
  1440. *
  1441. * FUNCTION
  1442. *
  1443. *   polysolve
  1444. *
  1445. * INPUT
  1446. *   
  1447. * OUTPUT
  1448. *   
  1449. * RETURNS
  1450. *   
  1451. * AUTHOR
  1452. *
  1453. *   Alexander Enzmann
  1454. *   
  1455. * DESCRIPTION
  1456. *
  1457. *   Root solver based on the Sturm sequences for a polynomial.
  1458. *
  1459. * CHANGES
  1460. *
  1461. *   Okt 1996 : Added bug fix by Heiko Eissfeldt. [DB]
  1462. *
  1463. ******************************************************************************/
  1464.  
  1465. static int polysolve(int order, DBL *Coeffs, DBL *roots)
  1466. {
  1467.   polynomial sseq[MAX_ORDER+1];
  1468.   DBL min_value, max_value;
  1469.   int i, nroots, np, atmin, atmax;
  1470.  
  1471.   /* Put the coefficients into the top of the stack. */
  1472.  
  1473.   for (i = 0; i <= order; i++)
  1474. {
  1475.     sseq[0].coef[order-i] = Coeffs[i] / Coeffs[0];
  1476.   }
  1477.  
  1478.   /* Build the Sturm sequence */
  1479.  
  1480.   np = buildsturm(order, &sseq[0]);
  1481.  
  1482.   /* Get the total number of visible roots */
  1483.  
  1484.   if ((nroots = visible_roots(np, sseq, &atmin, &atmax)) == 0)
  1485. {
  1486.     return(0);
  1487.   }
  1488.  
  1489.   /* Bracket the roots */
  1490.  
  1491.   min_value = 0.0;
  1492.   max_value = Max_Distance;
  1493.  
  1494.   atmin = numchanges(np, sseq, min_value);
  1495.   atmax = numchanges(np, sseq, max_value);
  1496.  
  1497.   nroots = atmin - atmax;
  1498.  
  1499.   if (nroots == 0)
  1500. {
  1501.     return(0);
  1502.   }
  1503.  
  1504.   /* perform the bisection. */
  1505.  
  1506.   return(sbisect(np, sseq, min_value, max_value, atmin, atmax, roots));
  1507. }
  1508.  
  1509.  
  1510.  
  1511. /*****************************************************************************
  1512. *
  1513. * FUNCTION
  1514. *
  1515. *   Solve_Polynomial
  1516. *
  1517. * INPUT
  1518. *
  1519. *   n       - order of polynomial
  1520. *   c       - coefficients
  1521. *   r       - roots
  1522. *   sturm   - TRUE, if sturm should be used for n=3,4
  1523. *   epsilon - Tolerance to discard small root
  1524. *   
  1525. * OUTPUT
  1526. *
  1527. *   r
  1528. *   
  1529. * RETURNS
  1530. *
  1531. *   int - number of roots found
  1532. *   
  1533. * AUTHOR
  1534. *
  1535. *   Dieter Bayer
  1536. *   
  1537. * DESCRIPTION
  1538. *
  1539. *   Solve the polynomial equation
  1540. *
  1541. *     c[0] * x ^ n + c[1] * x ^ (n-1) + ... + c[n-1] * x + c[n] = 0
  1542. *
  1543. *   If the equation has a root r, |r| < epsilon, this root is eliminated
  1544. *   and the equation of order n-1 will be solved. This will avoid the problem
  1545. *   of "surface acne" in (most) cases while increasing the speed of the
  1546. *   root solving (polynomial's order is reduced by one).
  1547. *
  1548. *   WARNING: This function can only be used for polynomials if small roots
  1549. *   (i.e. |x| < epsilon) are not needed. This is the case for ray/object
  1550. *   intersection tests because only intersecions with t > 0 are valid.
  1551. *
  1552. *   NOTE: Only one root at x = 0 will be eliminated.
  1553. *
  1554. *   NOTE: If epsilon = 0 no roots will be eliminated.
  1555. *
  1556. *
  1557. *   The method and idea for root elimination was taken from:
  1558. *
  1559. *     Han-Wen Nienhuys, "Polynomials", Ray Tracing News, July 6, 1994,
  1560. *     Volume 7, Number 3
  1561. *
  1562. *
  1563. * CHANGES
  1564. *
  1565. *   Jul 1994 : Creation.
  1566. *
  1567. ******************************************************************************/
  1568.  
  1569. int Solve_Polynomial(int n, DBL *c0, DBL *r, int sturm, DBL epsilon)
  1570. {
  1571.   int roots, i;
  1572.   DBL *c;
  1573.  
  1574.   Increase_Counter(stats[Polynomials_Tested]);
  1575.  
  1576.   roots = 0;
  1577.  
  1578.   /*
  1579.    * Determine the "real" order of the polynomial, i.e.
  1580.    * eliminate small leading coefficients.
  1581.    */
  1582.  
  1583.   i = 0;
  1584.  
  1585.   while ((fabs(c0[i]) < SMALL_ENOUGH) && (i < n))
  1586. {
  1587.     i++;
  1588.   }
  1589.  
  1590.   n -= i;
  1591.  
  1592.   c = &c0[i];
  1593.  
  1594.   switch (n)
  1595. {
  1596.     case 0:
  1597.  
  1598.       break;
  1599.  
  1600.     case 1:
  1601.  
  1602.       /* Solve linear polynomial. */
  1603.  
  1604.       if (c[0] != 0.0)
  1605.       {
  1606.         r[roots++] = -c[1] / c[0];
  1607.       }
  1608.  
  1609.       break;
  1610.  
  1611.     case 2:
  1612.  
  1613.       /* Solve quadratic polynomial. */
  1614.  
  1615.       roots = solve_quadratic(c, r);
  1616.  
  1617.       break;
  1618.  
  1619.     case 3:
  1620.  
  1621.       /* Root elimination? */
  1622.  
  1623.       if (epsilon > 0.0)
  1624.       {
  1625.         if ((c[2] != 0.0) && (fabs(c[3]/c[2]) < epsilon))
  1626.         {
  1627.           Increase_Counter(stats[Roots_Eliminated]);
  1628.  
  1629.           roots = solve_quadratic(c, r);
  1630.  
  1631.           break;
  1632.         }
  1633.       }
  1634.  
  1635.       /* Solve cubic polynomial. */
  1636.  
  1637.       if (sturm)
  1638.       {
  1639.         roots = polysolve(3, c, r);
  1640.       }
  1641.       else
  1642.       {
  1643.         roots = solve_cubic(c, r);
  1644.       }
  1645.  
  1646.       break;
  1647.  
  1648.     case 4:
  1649.  
  1650.       /* Root elimination? */
  1651.  
  1652.       if (epsilon > 0.0)
  1653.       {
  1654.         if ((c[3] != 0.0) && (fabs(c[4]/c[3]) < epsilon))
  1655.         {
  1656.           Increase_Counter(stats[Roots_Eliminated]);
  1657.  
  1658.           if (sturm)
  1659.           {
  1660.             roots = polysolve(3, c, r);
  1661.           }
  1662.           else
  1663.           {
  1664.             roots = solve_cubic(c, r);
  1665.           }
  1666.  
  1667.           break;
  1668.         }
  1669.       }
  1670.  
  1671.       /* Test for difficult coeffs. */
  1672.  
  1673.       if (difficult_coeffs(4, c))
  1674.       {
  1675.         sturm = TRUE;
  1676.       }
  1677.  
  1678.       /* Solve quartic polynomial. */
  1679.  
  1680.       if (sturm)
  1681.       {
  1682.         roots = polysolve(4, c, r);
  1683.       }
  1684.       else
  1685.       {
  1686.         roots = solve_quartic(c, r);
  1687.       }
  1688.  
  1689.       break;
  1690.  
  1691.     default:
  1692.  
  1693.       if (epsilon > 0.0)
  1694.       {
  1695.         if ((c[n-1] != 0.0) && (fabs(c[n]/c[n-1]) < epsilon))
  1696.         {
  1697.           Increase_Counter(stats[Roots_Eliminated]);
  1698.  
  1699.           roots = polysolve(n-1, c, r);
  1700.         }
  1701.       }
  1702.  
  1703.       /* Solve n-th order polynomial. */
  1704.  
  1705.       roots = polysolve(n, c, r);
  1706.  
  1707.       break;
  1708.   }
  1709.  
  1710.   return(roots);
  1711. }
  1712.