home *** CD-ROM | disk | FTP | other *** search
/ Stars of Shareware: Raytrace & Morphing / SOS-RAYTRACE.ISO / programm / source / rayce27s / solve.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-02-02  |  29.3 KB  |  1,360 lines
  1. /*
  2.  * solve.c -- solve equations
  3.  * 
  4.  * The larger part of this code (the sturm sequence code) is from Art, the
  5.  * raytracer from the rendering package Vort. It was written by David G
  6.  * Hook, Bernie Kirby and Peter McAree all the credits should go to them.
  7.  * You can reach them via <echidna@munnari.oz.au> The solve_quartic and
  8.  * solve_cubic are from Graphics Gems I "roots3and4.c", by Jochen
  9.  * Schwarze.
  10.  */
  11.  
  12. /*
  13.  * Ok...
  14.  * 
  15.  * It's not my code, but i'll try explain what this is all about.
  16.  * 
  17.  * This module solves polynomial equations, or equations that look like:
  18.  * 
  19.  * a_n x^n + a_(n-1) x^(n-1) + .. + a_0 = 0
  20.  * 
  21.  * n is called the order of the equation, is denoted by "ord" in the code.
  22.  * Note that any solution p(s) = 0 to the equation also is a root of the
  23.  * polynomial: if p(s) = 0 is a solution, then the polynomial can be
  24.  * written as p(x) = q(x) (x-s), with degree(q) = degree(p) - 1.
  25.  * 
  26.  * In this application we are only interested in so called "simple roots", ie
  27.  * if s is a root ( p(s) = 0 and p(x) = q(x) (x-s) ), then q(s) != 0.
  28.  * 
  29.  * This has  numerical reason: take (for example) p(x) := (x-1)^2. p clearly
  30.  * has a double root: p(1) = 0. But if we would start calculating,
  31.  * intermediary results would look like it came from p(x) = (x-1)^2 +
  32.  * 1e-20, which has no solutions. So double roots can be hard to detect.
  33.  * 
  34.  * If n <= 3, then an "exact" solution is possible. (Even with n == 4, an out
  35.  * of the box solution is possible. However, it is numerically unstable.
  36.  * has been proved that for n > 4, no general solution method exists.
  37.  * 
  38.  * If n >= 4, we have to resort to numerical techniques. The first attempt
  39.  * used was equation_solver(), which has been #undef'd ^H^H^Hdeleted,
  40.  * since it is 3x slower for quartics than sturm sequences.
  41.  * 
  42.  * It works like this: let us assume p(t) = 0 is the equation we want to
  43.  * solve.
  44.  * 
  45.  * 1. if degree <= 3, then solve exactly
  46.  * 
  47.  * 2. if degree > 3, then take derivative, (recursive) find it's roots with
  48.  * equation_solver().
  49.  * 
  50.  * These roots of the derivative are the extrema of p(t). We then calculate
  51.  * these extrema. Additionally, p(0) and p(LARGE) are boundary extrema.
  52.  * Any zero of p would lie between two extrema with different sign, and
  53.  * (vice versa) Between two extrema with different signs there always is a
  54.  * zero. We can find such a zero with binary chop.
  55.  * 
  56.  * The sturm sequence uses also this property: A sturm sequence can be used
  57.  * to detect sign changes of p(x) for x in the interval [a,b).
  58.  * 
  59.  * A sturm sequence is a sequence of polynomials,
  60.  * 
  61.  * p_0, p_1, ... p_m
  62.  * 
  63.  * because it's a sturm seq,  these polynomials have special properties. From
  64.  * these properties, it can be deduced that the number of sign changes in
  65.  * p_0(x) in the interval [a,b) is w(a) - w(b). Here, w(c) means the
  66.  * number of sign changes in
  67.  * 
  68.  * p_0(c), p_1(c), .. p_m(c)
  69.  * 
  70.  * A sturm sequence can easily be constructed using polynome division. For
  71.  * details: see J. Stoer, R. Bulirsch "Introduction to Numerical Analysis"
  72.  * pages 281 and further.
  73.  * 
  74.  * The number of sign changes between a and b equals the number of roots
  75.  * between a and b. So we could use binary chop too: after we've divided
  76.  * [a,b) into [a,c) and [c,b) we can check the number of signchanges in
  77.  * [a,c) and [c,b), until the number of signchanges is one, and then we're
  78.  * sure that one root is sitting in the interval.
  79.  */
  80.  
  81. #include "ray.h"
  82. #include "proto.h"
  83. #include "extern.h"
  84.  
  85. /*
  86.  * a coefficient smaller than SMALL_ENOUGH is considered to be zero (0.0).
  87.  */
  88.  
  89. /* epsilon surrounding for near zero values */
  90. #define     EQN_EPS     1e-9
  91.  
  92. /* cubic root. */
  93. #define     cbrt(x)     ((x) > 0.0 ? pow((double)(x), 1.0/3.0) : \
  94.                           ((x) < 0.0 ? -pow((double)-(x), 1.0/3.0) : 0.0))
  95.  
  96. #define    MAXP2        32
  97.  
  98. /* max number of iterations */
  99. #define          MAXIT             800
  100.  
  101. /* smallest relative error we want */
  102. #define    RELERROR    1.0e-10
  103.  
  104.  
  105. PRIVATE int     solve_cubic(double *coeffs, double *roots);
  106.  
  107. PRIVATE int     solve_rt_linear(double *coeff, double *roots);
  108. PRIVATE int     solve_rt_quartic(double *coef, double *roots);
  109. PRIVATE int     solve_rt_cubic(double *coeffs, double *roots);
  110. PRIVATE int     solve_rt_quadric(double *coeffs, double *roots);
  111.  
  112.  
  113. PRIVATE int     buildsturm(spoly *sseq);
  114. PRIVATE int     numchanges(int np, spoly *sseq, double a);
  115. PRIVATE int     numroots(int np, spoly *sseq, int *atneg, int *atpos);
  116. PRIVATE int     modp(spoly *u, spoly *v, spoly *r);
  117. PRIVATE int     modrf(int ord, double *coef, register double a, double b, double *val);
  118. PRIVATE int     evalatb(int np, register spoly *sseq, register double b);
  119. PRIVATE int     evalata(register int np, register spoly *sseq, register double a);
  120. PRIVATE int     evalat0(register int np, register spoly *sseq);
  121. PRIVATE int     sturm_solve(spoly *thepoly, double *roots, bool chkall);
  122. PRIVATE void    sbisect(int np, spoly *sseq, double min, double max, int atmin, int atmax, double *roots);
  123. PRIVATE int     csgsbisect(int np, spoly *sseq, double min, double max, int atmin, int atmax, double *roots);
  124.  
  125. /***************************************************************************
  126.  * ENTRY POINT: solve polynomial, return number of nonnegative roots       *
  127.  ***************************************************************************/
  128.  
  129. PUBLIC int
  130. solve_rt_poly(spoly *thepoly, double *roots, bool chkall)
  131. {
  132.     int             r;
  133.  
  134.     /* avoid degenerate cases. */
  135.     poly_clean(thepoly);
  136.  
  137.  
  138.     /* special cases can be handled much faster: */
  139.     switch (thepoly->ord) {
  140.     case 4:
  141.     return solve_rt_quartic(thepoly->coef, roots);
  142.     case 2:
  143.     return solve_rt_quadric(thepoly->coef, roots);
  144.     case 3:
  145.     return solve_rt_cubic(thepoly->coef, roots);
  146.     case 1:
  147.     return solve_rt_linear(thepoly->coef, roots);
  148.     case 0:
  149.     return 0;
  150.     }
  151.  
  152.     r = sturm_solve(thepoly, roots, chkall);
  153.  
  154.     return r;
  155. }
  156.  
  157. PUBLIC int
  158. sturm_solve_rt_poly(spoly *thepoly, double *roots, bool chkall)
  159. {
  160.     int             r;
  161.  
  162.     /* not really necessary for these polys. */
  163.     if (thepoly->ord <= 2)
  164.     return solve_rt_poly(thepoly, roots, chkall);
  165.  
  166.     /* avoid degenerate cases. */
  167.     poly_clean(thepoly);
  168.  
  169.     while (thepoly->ord > 0 && thepoly->coef[thepoly->ord] == 0.0)
  170.     thepoly->ord--;
  171.  
  172.     r = sturm_solve(thepoly, roots, chkall);
  173.     return r;
  174. }
  175.  
  176. /*
  177.  * Polysolver, based on Eric's Sturm sequence polysolver (function below
  178.  * was taken from torus.c
  179.  */
  180. PRIVATE int
  181. sturm_solve(spoly *thepoly, double *roots, bool chkall)
  182. {
  183.     int             np,        /* number of polys in the sturm sequence */
  184.                     numroots,
  185.                     nroots,
  186.                     atmin,
  187.                     atinf;
  188.  
  189.     spoly           sseq[MAX_ORDER];    /* structure to hold the sequence */
  190.  
  191.     register int    its,
  192.                     i;
  193.     register float  min,
  194.                     max;
  195.  
  196.     /* copy polynomial */
  197.     for (i = 0; i <= thepoly->ord; i++) {
  198.     sseq[0].coef[i] = thepoly->coef[i];
  199.     }
  200.  
  201.     sseq[0].ord = thepoly->ord;
  202.  
  203.     /* build it's sturm sequence. */
  204.     np = buildsturm(sseq);
  205.  
  206.     /* find sign changes in the sequence at 0 */
  207.     atmin = evalat0(np, sseq);
  208.  
  209.     /* find sign changes in the sequence at infinity */
  210.     atinf = evalatb(np, sseq, 0.0);
  211.  
  212.     /*
  213.      * the number of sign changes of thepoly in [a,b] is atmin - atinf.
  214.      * This is the important thing of Sturm seqs
  215.      */
  216.     nroots = atmin - atinf;
  217.  
  218.     /* no roots which are bigger than 0 */
  219.     if (nroots == 0) {
  220.     return 0;
  221.     }
  222.     /* hmm ... I wonder why. */
  223.     atmin = evalata(np, sseq, tolerance);
  224.  
  225.     nroots = atmin - atinf;
  226.  
  227.     if (nroots == 0) {
  228.     return 0;
  229.     }
  230.     min = tolerance;
  231.     max = min + 2.0;
  232.  
  233.     /*
  234.      * find the range where all those signchanges occur, by doubling the
  235.      * interval size.
  236.      */
  237.     for (its = 0; its < MAXP2; its++) {
  238.     atinf = evalatb(np, sseq, max);
  239.  
  240.     numroots = atmin - atinf;
  241.  
  242.     if (numroots == nroots)    /* [min, max] covers all signchanges,
  243.                  * therefore it covers all zero's */
  244.         break;
  245.  
  246.     max *= 2.0;
  247.     }
  248.  
  249.  
  250.     if (nroots == 1 || !chkall) {
  251.     /* only one root, (or we're only interested in the smallest) */
  252.     sbisect(np, sseq, min, max, atmin, atinf, roots);
  253.     return 1;
  254.  
  255.     } else {            /* we'll need this if CSG is implemented. */
  256.     if (csgsbisect(np, sseq, min, max, atmin, atinf, roots) != nroots)
  257.         assert(FALSE);
  258.     return nroots;
  259.     }
  260. }
  261.  
  262.  
  263. /*
  264.  * evalat0
  265.  * 
  266.  * calculate the sign changes in the sturm sequence in smat at 0
  267.  */
  268. PRIVATE int
  269. evalat0(register int np, register spoly *sseq)
  270. {
  271.     register int    at0;
  272.     register spoly *s,
  273.                    *ls;
  274.  
  275.     at0 = 0;
  276.  
  277.     ls = sseq;
  278.     for (s = sseq + 1; s <= sseq + np; s++) {
  279.     if (ls->coef[0] * s->coef[0] < 0 || (ls->coef[0] == 0.0 && s->coef[0] != 0.0))
  280.         at0++;
  281.     ls = s;
  282.     }
  283.  
  284.     return (at0);
  285. }
  286.  
  287. /*
  288.  * evalata
  289.  * 
  290.  * calculate the sign changes in the sturm sequence in smat at min value a
  291.  */
  292. PRIVATE int
  293. evalata(register int np, register spoly *sseq, register double a)
  294. {
  295.     register int    ata;
  296.     register double f,
  297.                     lf;
  298.     register spoly *s;
  299.  
  300.     ata = 0;
  301.  
  302.     lf = poly_eval(a, sseq[0].ord, sseq[0].coef);
  303.  
  304.     for (s = sseq + 1; s <= sseq + np; s++) {
  305.     f = poly_eval(a, s->ord, s->coef);
  306.     if (lf * f < 0 || lf == 0.0)
  307.         ata++;
  308.     lf = f;
  309.     }
  310.  
  311.     return (ata);
  312. }
  313.  
  314. /*
  315.  * evalatb
  316.  * 
  317.  * calculate the sign changes in the sturm sequence in smat at max value b
  318.  */
  319. PRIVATE int
  320. evalatb(int np, register spoly *sseq, register double b)
  321. {
  322.     int             atb;
  323.     register double f,
  324.                     lf;
  325.     register spoly *s,
  326.                    *ls;
  327.  
  328.     atb = 0;
  329.  
  330.     if (b == 0.0) {
  331.     ls = sseq;
  332.     for (s = sseq + 1; s <= sseq + np; s++) {
  333.         if (ls->coef[ls->ord] * s->coef[s->ord] < 0 || (ls->coef[ls->ord] == 0.0 && s->coef[s->ord] != 0.0))
  334.         atb++;
  335.         ls = s;
  336.     }
  337.     } else {
  338.     lf = poly_eval(b, sseq[0].ord, sseq[0].coef);
  339.     for (s = sseq + 1; s <= sseq + np; s++) {
  340.         f = poly_eval(b, s->ord, s->coef);
  341.         if (lf * f < 0 || lf == 0.0)
  342.         atb++;
  343.         lf = f;
  344.     }
  345.     }
  346.  
  347.     return (atb);
  348. }
  349.  
  350.  
  351. /*
  352.  * modrf
  353.  * 
  354.  * uses the modified regula-falsi method to evaluate the root in interval
  355.  * [a,b] of the polynomial described in coef. The root is returned is
  356.  * returned in *val. The routine returns zero if it can't converge.
  357.  */
  358. PRIVATE int
  359. modrf(int ord, double *coef, register double a, double b, double *val)
  360. {
  361.     register int    its;
  362.     register double fa,
  363.                     fb,
  364.                     x,
  365.                     fx,
  366.                     lfx;
  367.     register double *fp,
  368.                    *scoef,
  369.                    *ecoef;
  370.  
  371.     scoef = coef;
  372.     ecoef = &coef[ord];
  373.  
  374.     fb = fa = *ecoef;
  375.     for (fp = ecoef - 1; fp >= scoef; fp--) {
  376.     fa = a * fa + *fp;
  377.     fb = b * fb + *fp;
  378.     }
  379.  
  380.     /*
  381.      * if there is no sign difference the method won't work
  382.      */
  383.     if (fa * fb > 0.0)
  384.     return (0);
  385.  
  386.     if (fabs(fa) < RELERROR) {
  387.     *val = a;
  388.     return (1);
  389.     }
  390.     if (fabs(fb) < RELERROR) {
  391.     *val = b;
  392.     return (1);
  393.     }
  394.     lfx = fa;
  395.  
  396.     for (its = 0; its < MAXIT; its++) {
  397.  
  398.     x = (fb * a - fa * b) / (fb - fa);
  399.  
  400.     fx = *ecoef;
  401.     for (fp = ecoef - 1; fp >= scoef; fp--)
  402.         fx = x * fx + *fp;
  403.  
  404.     if (fabs(x) > RELERROR) {
  405.         if (fabs(fx / x) < RELERROR) {
  406.         *val = x;
  407.         return (1);
  408.         }
  409.     } else if (fabs(fx) < RELERROR) {
  410.         *val = x;
  411.         return (1);
  412.     }
  413.     if ((fa * fx) < 0) {
  414.         b = x;
  415.         fb = fx;
  416.         if ((lfx * fx) > 0)
  417.         fa /= 2;
  418.     } else {
  419.         a = x;
  420.         fa = fx;
  421.         if ((lfx * fx) > 0)
  422.         fb /= 2;
  423.     }
  424.  
  425.     lfx = fx;
  426.     }
  427.  
  428.     if (fabs(fx) > ALG_TOLERANCE) {
  429. #ifdef DEBUG
  430.     if (debug_options & DEBUGRUNTIME)
  431.         warning("modrf overflow %f %f %f\n", a, b, fx);
  432. #endif
  433.     }
  434.     *val = x;
  435.  
  436.     return (1);
  437. }
  438.  
  439. /*
  440.  * sbisect
  441.  * 
  442.  * uses a bisection based on the sturm sequence for the polynomial described
  443.  * in sseq to isolate the lowest positve root of the polynomial in the
  444.  * bounds min and max. The root is returned in *roots.
  445.  */
  446. PRIVATE void
  447. sbisect(int np, spoly *sseq, double min, double max, int atmin, int atmax, double *roots)
  448. {
  449.     register double mid;
  450.     register int    n1,
  451.                     its,
  452.                     atmid;
  453.  
  454.  
  455.     for (its = 0; its < MAXIT; its++) {
  456.  
  457.     mid = (min + max) / 2;
  458.  
  459.     atmid = numchanges(np, sseq, mid);
  460.  
  461.     n1 = atmin - atmid;
  462.  
  463.     if (n1 == 1) {
  464.         /*
  465.          * first try a less expensive technique.
  466.          */
  467.         if (modrf(sseq->ord, sseq->coef, min, mid, &roots[0]))
  468.         return;
  469.  
  470.         /*
  471.          * if we get here we have to evaluate the root the hard way by
  472.          * using the Sturm sequence. It also means the root is
  473.          * repeated an even number of times.
  474.          */
  475.         for (its = 0; its < MAXIT; its++) {
  476.         mid = (min + max) / 2;
  477.  
  478.         atmid = numchanges(np, sseq, mid);
  479.  
  480.         if (fabs(mid) > RELERROR) {
  481.             if (fabs((max - min) / mid) < RELERROR) {
  482.             roots[0] = mid;
  483.             return;
  484.             }
  485.         } else if (fabs(max - min) < RELERROR) {
  486.             roots[0] = mid;
  487.             return;
  488.         }
  489.         if ((atmin - atmid) == 0)
  490.             min = mid;
  491.         else
  492.             max = mid;
  493.         }
  494.  
  495. #ifdef DEBUG
  496.         if (its == MAXIT && (max - min) > ALG_TOLERANCE &&
  497.         (debug_options & DEBUGRUNTIME)) {
  498.         warning("sbisect: overflow min %f max %f diff %e\n", min, max, max - min);
  499.         }
  500. #endif
  501.         roots[0] = mid;
  502.  
  503.         return;
  504.     } else if (n1 == 0)
  505.         min = mid;
  506.     else
  507.         max = mid;
  508.  
  509.     }
  510.  
  511.     roots[0] = min;        /* multiple roots close together */
  512. #ifdef DEBUG
  513.     if (its == MAXIT && (max - min) > ALG_TOLERANCE && (debug_options & DEBUGRUNTIME))
  514.     warning("sbisect: overflow min %f max %f diff %e\n", min, max, max - min);
  515. #endif
  516. }
  517.  
  518. /*
  519.  * csgsbisect
  520.  * 
  521.  * uses a bisection based on the sturm sequence for the polynomial described
  522.  * in sseq to find the positive roots in the bounds min to max, the roots
  523.  * are returned in the roots array in order of magnitude.
  524.  */
  525. PRIVATE int
  526. csgsbisect(int np, spoly *sseq, double min, double max, int atmin, int atmax, double *roots)
  527. {
  528.     double          mid;
  529.  
  530.     int             n1,
  531.                     n2,
  532.                     its,
  533.                     atmid,
  534.                     nroot,
  535.                     nr;
  536.  
  537.     /* assume the worst */
  538.     if ((nroot = atmin - atmax) == 1) {
  539.  
  540.     /*
  541.      * first try a less expensive technique.
  542.      */
  543.     if (modrf(sseq->ord, sseq->coef, min, max, &roots[0])) {
  544.         return 1;
  545.     }
  546.     /*
  547.      * if we get here we have to evaluate the root the hard way by
  548.      * using the Sturm sequence. It also means the root is repeated an
  549.      * even number of times.
  550.      */
  551.     for (its = 0; its < MAXIT; its++) {
  552.         mid = (min + max) / 2;
  553.  
  554.         atmid = numchanges(np, sseq, mid);
  555.  
  556.         if (fabs(mid) > RELERROR) {
  557.         if (fabs((max - min) / mid) < RELERROR)
  558.             break;
  559.         } else if (fabs(max - min) < RELERROR)
  560.         break;
  561.  
  562.         if ((atmin - atmid) == 0)
  563.         min = mid;
  564.         else
  565.         max = mid;
  566.     }
  567.  
  568. #ifdef DEBUG
  569.     if (its == MAXIT && (debug_options & DEBUGRUNTIME))
  570.         warning("csgsbisect: overflow min %f max %f diff %e\n", min, max, max - min);
  571.  
  572. #endif
  573.     roots[1] = roots[0] = mid;
  574.     return 2;
  575.     }
  576.     /*
  577.      * more than one root in the interval, we have to bisect...
  578.      */
  579.     nr = 0;
  580.     for (its = 0; its < MAXIT; its++) {
  581.  
  582.     mid = (min + max) / 2;
  583.  
  584.     atmid = numchanges(np, sseq, mid);
  585.  
  586.     n1 = atmin - atmid;
  587.     n2 = atmid - atmax;
  588.  
  589.     if (n1 != 0 && n2 != 0) {
  590.         nr += csgsbisect(np, sseq, min, mid, atmin, atmid, roots);
  591.         nr += csgsbisect(np, sseq, mid, max, atmid, atmax, &roots[n1]);
  592.         break;
  593.     }
  594.     if (n1 == 0)
  595.         min = mid;
  596.     else
  597.         max = mid;
  598.     }
  599.  
  600.     if (its == MAXIT) {
  601.     int             i;
  602.  
  603. #ifdef DEBUG
  604.     if (debug_options & DEBUGRUNTIME) {
  605.         warning("csgsbisect: roots too close together\n");
  606.         warning("csgsbisect: overflow min %f max %f diff %e nroot %d n1 %d n2 %d\n", min, max, max - min, nroot, n1, n2);
  607.     }
  608. #endif
  609.  
  610.     n1 = atmin - atmax;
  611.     nr += n1;
  612.     for (i = 0; i < n1; i++)
  613.         roots[i] = mid;
  614.     }
  615.     return nr;
  616. }
  617.  
  618. /*
  619.  * modp
  620.  * 
  621.  * calculates the modulus of u(x) / v(x) leaving it in r, it returns 0 if
  622.  * r(x) is a constant.
  623.  * 
  624.  * note: this function assumes the leading coefficient of v is 1 or -1
  625.  */
  626. PRIVATE int
  627. modp(spoly *u, spoly *v, spoly *r)
  628. {
  629.     int             k,
  630.                     j;
  631.     double         *nr,
  632.                    *end,
  633.                    *uc;
  634.  
  635.     nr = r->coef;
  636.     end = &u->coef[u->ord];
  637.  
  638.     uc = u->coef;
  639.     while (uc <= end)
  640.     *nr++ = *uc++;
  641.  
  642.     if (v->coef[v->ord] < 0.0) {
  643.  
  644.     for (k = u->ord - v->ord - 1; k >= 0; k -= 2)
  645.         r->coef[k] = -r->coef[k];
  646.  
  647.     for (k = u->ord - v->ord; k >= 0; k--)
  648.         for (j = v->ord + k - 1; j >= k; j--)
  649.         r->coef[j] = -r->coef[j] - r->coef[v->ord + k] * v->coef[j - k];
  650.     } else {
  651.     for (k = u->ord - v->ord; k >= 0; k--)
  652.         for (j = v->ord + k - 1; j >= k; j--)
  653.         r->coef[j] -= r->coef[v->ord + k] * v->coef[j - k];
  654.     }
  655.  
  656.     k = v->ord - 1;
  657.     while (k >= 0 && r->coef[k] == 0.0)
  658.     k--;
  659.  
  660.     r->ord = (k < 0) ? 0 : k;
  661.  
  662.     return (r->ord);
  663. }
  664.  
  665. /*
  666.  * buildsturm
  667.  * 
  668.  * build up a sturm sequence for a polynomial in smat, returning the number
  669.  * of polynomials in the sequence
  670.  */
  671. PRIVATE int
  672. buildsturm(spoly *sseq)
  673. {
  674.     int             i,
  675.                     ord;
  676.     double          f,
  677.                    *fp,
  678.                    *fc;
  679.     spoly          *sp;
  680.  
  681.     ord = sseq[0].ord;
  682.     sseq[1].ord = ord - 1;
  683.  
  684.     /*
  685.      * calculate the derivative and normalise the leading coefficient.
  686.      */
  687.     f = fabs(sseq[0].coef[ord] * ord);
  688.     fp = sseq[1].coef;
  689.     fc = sseq[0].coef + 1;
  690.     for (i = 1; i <= ord; i++)
  691.     *fp++ = *fc++ * i / f;
  692.  
  693.     /*
  694.      * construct the rest of the Sturm sequence
  695.      */
  696.     for (sp = sseq + 2; modp(sp - 2, sp - 1, sp); sp++) {
  697.  
  698.     /*
  699.      * reverse the sign and normalise
  700.      */
  701.     f = -fabs(sp->coef[sp->ord]);
  702.     for (fp = &sp->coef[sp->ord]; fp >= sp->coef; fp--)
  703.         *fp /= f;
  704.     }
  705.  
  706.     sp->coef[0] = -sp->coef[0];    /* reverse the sign */
  707.  
  708.     return (sp - sseq);
  709. }
  710.  
  711. /*
  712.  * numroots
  713.  * 
  714.  * return the number of distinct real roots of the polynomial described in
  715.  * sseq.
  716.  */
  717. PRIVATE int
  718. numroots(int np,
  719.      spoly *sseq,
  720.      int *atneg, int *atpos)
  721. {
  722.     int             atposinf,
  723.                     atneginf;
  724.     spoly          *s;
  725.     double          f,
  726.                     lf;
  727.  
  728.     atposinf = atneginf = 0;
  729.  
  730.     /*
  731.      * changes at positve infinity
  732.      */
  733.     lf = sseq[0].coef[sseq[0].ord];
  734.  
  735.     for (s = sseq + 1; s <= sseq + np; s++) {
  736.     f = s->coef[s->ord];
  737.     if (lf == 0.0 || lf * f < 0)
  738.         atposinf++;
  739.     lf = f;
  740.     }
  741.  
  742.     /*
  743.      * changes at negative infinity
  744.      */
  745.     if (sseq[0].ord & 1)
  746.     lf = -sseq[0].coef[sseq[0].ord];
  747.     else
  748.     lf = sseq[0].coef[sseq[0].ord];
  749.  
  750.     for (s = sseq + 1; s <= sseq + np; s++) {
  751.     if (s->ord & 1)
  752.         f = -s->coef[s->ord];
  753.     else
  754.         f = s->coef[s->ord];
  755.     if (lf == 0.0 || lf * f < 0)
  756.         atneginf++;
  757.     lf = f;
  758.     }
  759.  
  760.     *atneg = atneginf;
  761.     *atpos = atposinf;
  762.  
  763.     return (atneginf - atposinf);
  764. }
  765.  
  766. /*
  767.  * numchanges
  768.  * 
  769.  * return the number of sign changes in the Sturm sequence in sseq at the
  770.  * value a.
  771.  */
  772. PRIVATE int
  773. numchanges(int np,
  774.        spoly *sseq,
  775.        double a)
  776. {
  777.     int             changes;
  778.     double          f,
  779.                     lf;
  780.     spoly          *s;
  781.  
  782.     changes = 0;
  783.  
  784.     lf = poly_eval(a, sseq[0].ord, sseq[0].coef);
  785.  
  786.     for (s = sseq + 1; s <= sseq + np; s++) {
  787.     f = poly_eval(a, s->ord, s->coef);
  788.     if (lf == 0.0 || lf * f < 0)
  789.         changes++;
  790.     lf = f;
  791.     }
  792.  
  793.     return (changes);
  794. }
  795.  
  796.  
  797. /* solve linear equation */
  798. PRIVATE int
  799. solve_rt_linear(double *coeff, double *roots)
  800. {
  801.     if ((*roots = -coeff[0] / coeff[1]) > tolerance)
  802.     return 1;
  803.     else
  804.     return 0;
  805. }
  806.  
  807. /*
  808.  * SPECIAL ROOTFINDERS: quadric, cubic and quartic.
  809.  * 
  810.  * written by Jochen Schwarze and me from Graphics Gems I, by Jochen
  811.  * Schwarze, (schwarze@isa.de) quadratic: both by Jochen Schwarze,  from
  812.  * graphics gems cubic: one from me, and one by Jochen Schwarze quartic:
  813.  * by Jochen Schwarze.
  814.  * 
  815.  * 
  816.  * Here is the original comment header:
  817.  * 
  818.  * Roots3And4.c
  819.  * 
  820.  * Utility functions to find cubic and quartic roots, coefficients are passed
  821.  * like this:
  822.  * 
  823.  * c[0] + c[1]*x + c[2]*x^2 + c[3]*x^3 + c[4]*x^4 = 0
  824.  * 
  825.  * The functions return the number of non-complex roots and put the values
  826.  * into the s array.
  827.  * 
  828.  * Author:         Jochen Schwarze (schwarze@isa.de)
  829.  * 
  830.  * Jan 26, 1990    Version for Graphics Gems Oct 11, 1990    Fixed sign
  831.  * problem for negative q's in SolveQuartic (reported by Mark Podlipec),
  832.  * Old-style function definitions, IsZero() as a macro Nov 23, 1990 Some
  833.  * systems do not declare acos() and cbrt() in <math.h>, though the
  834.  * functions exist in the library. If large coefficients are used, EQN_EPS
  835.  * should be reduced considerably (e.g. to 1E-30), results will be correct
  836.  * but multiple roots might be reported more than once.
  837.  */
  838.  
  839.  
  840.  
  841.  
  842. /*
  843.  * all roots of quadratic eqn.
  844.  */
  845. int
  846. solve_quadric(double c[3], double s[2])
  847. {
  848.     double          p,
  849.                     q,
  850.                     D;
  851.  
  852.     /* normal form: x^2 + px + q = 0 */
  853.  
  854.     p = c[1] / (2 * c[2]);
  855.     q = c[0] / c[2];
  856.  
  857.     D = p * p - q;
  858.  
  859.     if (D < 0) {
  860.     return 0;
  861.     } else {
  862.     D = sqrt(D);
  863.  
  864.     s[0] = D - p;
  865.     s[1] = -D - p;
  866.     return 2;
  867.     }
  868. }
  869.  
  870. /*
  871.  * solve quadric, return roots sorted and bigger than 0 double roots are
  872.  * returned double, so CSG still works.
  873.  */
  874. int
  875. solve_rt_quadric(double *coeffs, double *roots)
  876. {
  877.     double          p,
  878.                     q,
  879.                     D,
  880.                     r1,
  881.                     r2;
  882.  
  883.  
  884.     /*
  885.      * pre: coeffs[2] != 0
  886.      */
  887.  
  888.     p = coeffs[1] / (2 * coeffs[2]);
  889.     q = coeffs[0] / coeffs[2];
  890.  
  891.     D = p * p - q;
  892.  
  893.     if (D < 0) {
  894.     return 0;
  895.     }
  896.     /*
  897.      * don't check for double roots. If there are any, we still need both
  898.      * of them.
  899.      */
  900.  
  901.     D = sqrt(D);
  902.     r1 = D - p;
  903.     r2 = -D - p;
  904.  
  905.     if (r1 > r2) {
  906.     D = r2;
  907.     r2 = r1;
  908.     r1 = D;
  909.     }
  910.     /* r2 is maximum. If r2 < epsilon, then shake it */
  911.     if (r2 < tolerance)
  912.     return 0;
  913.  
  914.     /* r2 > tolerance, but is r1 > tolerance? */
  915.     if (r1 < tolerance) {
  916.     *roots = r2;
  917.     return 1;
  918.     } else {
  919.     *roots++ = r1;
  920.     *roots = r2;
  921.     return 2;
  922.     }
  923.  
  924. }
  925.  
  926. /*
  927.  * solve an equation using Cardano's method. Numerically unstable with
  928.  * double roots.
  929.  * 
  930.  * How does it work?
  931.  * 
  932.  * Suppose we have an cubic equation
  933.  * 
  934.  * ax^3+ bx^2 + cx + d = 0
  935.  * 
  936.  * After division by a, we get
  937.  * 
  938.  * x^3 + b' x^2 + c' x + d' = 0
  939.  * 
  940.  * After substitution x = y -b'/3. We get
  941.  * 
  942.  * y^3 = 3py - 2q        (1)
  943.  * 
  944.  * For certain p and q. This is the big Cardano trick: we could solve the
  945.  * equation if we would find u and v such that
  946.  * 
  947.  * u^3 + v^3  = -2q and uv = p            (2)
  948.  * 
  949.  * since (u+v)^3 = u^3 + 3u^2v + 3v^2u + v^3. = 3(u+v) uv + u^3 + v^3
  950.  * 
  951.  * and y would be y = u+v. (2) turns out to be a quadratic equation, and it
  952.  * gives us one solution for (1).  If (2) has no solution (discriminant
  953.  * 4q^2 - 4p^3 < 0), an other method has to used. We can safely assume
  954.  * that y = 2sqrt(p) cos(theta), and (1) turns out to be
  955.  * 
  956.  * cos(3 theta) = -q/(p*(sqrt(p)))
  957.  * 
  958.  * This gives us 3 solutions.
  959.  */
  960.  
  961.  
  962. /*
  963.  * this one is for RT applications, it returns number of roots bigger than
  964.  * EPSILON, and the roots itself sorted in an array.
  965.  */
  966.  
  967. #ifdef UNDEFINED
  968. PRIVATE int
  969. solve_rt_cubic(double *coeffs, double *roots)
  970. {
  971.     double          tcoeffs[4],
  972.                     a,
  973.                     lin,
  974.                     cons;
  975.     int             noroots,
  976.                     i;
  977.  
  978.     int             n,
  979.                     j;
  980.     double          r[4];
  981.  
  982.     n = solve_quartic(coeffs, r);
  983.     qsort(r, n, sizeof(double), dcomp);
  984.  
  985.     j = 0;
  986.     for (i = 0; i < n; i++)
  987.     if (r[i] > tolerance)
  988.         roots[j++] = r[i];
  989.     return j;
  990.  
  991.  
  992.  
  993.     a = coeffs[3];
  994.     for (i = 0; i < 4; i++)
  995.     tcoeffs[i] = coeffs[i] / a;
  996.  
  997.     a = tcoeffs[2];
  998.  
  999.     lin = -sqr(a) / 3.0 + tcoeffs[1];
  1000.     cons = (a * a * a * 2) / 27.0 - a * tcoeffs[1] / 3.0 + tcoeffs[0];
  1001.  
  1002.     {
  1003.     double          p,
  1004.                     q,
  1005.                     u,
  1006.                     v,
  1007.                     d;
  1008.  
  1009.     p = -lin / 3;
  1010.     q = cons / 2;
  1011.  
  1012.     d = 4 * q * q - 4 * p * p * p;
  1013.  
  1014.     if (d >= 0) {
  1015.         d = sqrt(d);
  1016.         u = (-2 * q + d) / 2;
  1017.         v = (-2 * q - d) / 2;
  1018.  
  1019.         u = cbrt(u);
  1020.         v = cbrt(v);
  1021.  
  1022.         *roots = u + v - a / 3;
  1023.  
  1024.         if (*roots > tolerance)
  1025.         return 1;
  1026.         else
  1027.         return 0;
  1028.     } else {        /* discr < 0, use Vi\`ete's method */
  1029.         double          r,
  1030.                        *rp,
  1031.                         rt[3];
  1032.  
  1033.         noroots = 0;
  1034.         rp = rt;
  1035.  
  1036.         a = a / 3;
  1037.  
  1038.         d = acos(-q / (p * sqrt(p))) / 3;
  1039.         p = 2 * sqrt(p);
  1040.  
  1041.         if ((r = p * cos(d) - a) > tolerance)
  1042.         *rp++ = r;
  1043.  
  1044.         d += 2 * M_PI / 3;
  1045.         if ((r = p * cos(d) - a) > tolerance)
  1046.         *rp++ = r;
  1047.  
  1048.         d += 2 * M_PI / 3;
  1049.  
  1050.         if ((r = p * cos(d) - a) > tolerance)
  1051.         *rp++ = r;
  1052.  
  1053.  
  1054.         /*
  1055.          * calculation is done. Sort the roots. Smallest must be
  1056.          * first.
  1057.          */
  1058.  
  1059.         noroots = rp - rt;
  1060.  
  1061.         if (noroots <= 1) {
  1062.         *roots = rt[0];
  1063.         return noroots;
  1064.         }
  1065.         if (rt[0] > rt[1]) {
  1066.         r = rt[0];
  1067.         rt[0] = rt[1];
  1068.         rt[1] = r;
  1069.         }
  1070.         if (noroots == 2) {
  1071.         *roots++ = rt[0];
  1072.         *roots++ = rt[1];
  1073.         return 2;
  1074.         }
  1075.         if (rt[0] > rt[2]) {
  1076.         r = rt[0];
  1077.         rt[0] = rt[2];
  1078.         rt[2] = r;
  1079.         }
  1080.         /* now, rt[0] is the smallest */
  1081.         if (rt[1] > rt[2]) {
  1082.         r = rt[1];
  1083.         rt[1] = rt[2];
  1084.         rt[2] = r;
  1085.         }
  1086.         /* copy the roots */
  1087.         *roots++ = rt[0];
  1088.         *roots++ = rt[1];
  1089.         *roots++ = rt[2];
  1090.  
  1091.         return 3;
  1092.     }
  1093.     }
  1094. }
  1095.  
  1096. #endif
  1097.  
  1098. /*
  1099.  * solve cubic, same method as above. Don't sort roots; return all roots.
  1100.  */
  1101. PRIVATE int
  1102. solve_cubic(double c[4], double s[3])
  1103. {
  1104.     int             i,
  1105.                     num;
  1106.     double          sub;
  1107.     double          A,
  1108.                     B,
  1109.                     C;
  1110.     double          sq_A,
  1111.                     p,
  1112.                     q;
  1113.     double          cb_p,
  1114.                     D;
  1115.  
  1116.     /* normal form: x^3 + Ax^2 + Bx + C = 0 */
  1117.  
  1118.     A = c[2] / c[3];
  1119.     B = c[1] / c[3];
  1120.     C = c[0] / c[3];
  1121.  
  1122.     /*
  1123.      * substitute x = y - A/3 to eliminate quadric term: x^3 +px + q = 0
  1124.      */
  1125.  
  1126.     sq_A = A * A;
  1127.     p = 1.0 / 3 * (-1.0 / 3 * sq_A + B);
  1128.     q = 1.0 / 2 * (2.0 / 27 * A * sq_A - 1.0 / 3 * A * B + C);
  1129.  
  1130.     /* use Cardano's formula */
  1131.  
  1132.     cb_p = p * p * p;
  1133.     D = q * q + cb_p;
  1134.  
  1135.     if (ISZERO(D)) {
  1136.     if (ISZERO(q)) {    /* one triple solution */
  1137.         s[0] = 0;
  1138.         num = 1;
  1139.     } else {        /* one single and one double solution */
  1140.         double          u = cbrt(-q);
  1141.  
  1142.         s[0] = 2 * u;
  1143.         s[1] = -u;
  1144.         num = 2;
  1145.     }
  1146.     } else if (D < 0) {        /* Casus irreducibilis: three real
  1147.                  * solutions */
  1148.     double          phi = 1.0 / 3 * acos(-q / sqrt(-cb_p));
  1149.     double          t = 2 * sqrt(-p);
  1150.  
  1151.     s[0] = t * cos(phi);
  1152.     s[1] = -t * cos(phi + M_PI / 3);
  1153.     s[2] = -t * cos(phi - M_PI / 3);
  1154.     num = 3;
  1155.     } else {            /* one real solution */
  1156.     double          sqrt_D = sqrt(D);
  1157.     double          u = cbrt(sqrt_D - q);
  1158.     double          v = -cbrt(sqrt_D + q);
  1159.  
  1160.     s[0] = u + v;
  1161.     num = 1;
  1162.     }
  1163.  
  1164.     /* resubstitute */
  1165.  
  1166.     sub = 1.0 / 3 * A;
  1167.  
  1168.     for (i = 0; i < num; ++i)
  1169.     s[i] -= sub;
  1170.  
  1171.     return num;
  1172. }
  1173.  
  1174. /*
  1175.  * solve a quartic eqn. Ferrari was the first to think of this. (I believe
  1176.  * or was it Descartes?
  1177.  */
  1178.  
  1179. PRIVATE int
  1180. solve_quartic(double c[5], double s[4])
  1181. {
  1182.     double          coeffs[4];
  1183.     double          z,
  1184.                     u,
  1185.                     v,
  1186.                     sub;
  1187.     double          A,
  1188.                     B,
  1189.                     C,
  1190.                     D;
  1191.     double          sq_A,
  1192.                     p,
  1193.                     q,
  1194.                     r;
  1195.     int             i,
  1196.                     num;
  1197.  
  1198.     /* normal form: x^4 + Ax^3 + Bx^2 + Cx + D = 0 */
  1199.  
  1200.     A = c[3] / c[4];
  1201.     B = c[2] / c[4];
  1202.     C = c[1] / c[4];
  1203.     D = c[0] / c[4];
  1204.  
  1205.     /*
  1206.      * substitute x = y - A/4 to eliminate cubic term: x^4 + px^2 + qx + r
  1207.      * = 0
  1208.      */
  1209.  
  1210.     sq_A = A * A;
  1211.     p = -3.0 / 8 * sq_A + B;
  1212.     q = 1.0 / 8 * sq_A * A - 1.0 / 2 * A * B + C;
  1213.     r = -3.0 / 256 * sq_A * sq_A + 1.0 / 16 * sq_A * B - 1.0 / 4 * A * C + D;
  1214.  
  1215.     if (ISZERO(r)) {
  1216.     /* no absolute term: y(y^3 + py + q) = 0 */
  1217.  
  1218.     coeffs[0] = q;
  1219.     coeffs[1] = p;
  1220.     coeffs[2] = 0;
  1221.     coeffs[3] = 1;
  1222.  
  1223.     num = solve_cubic(coeffs, s);
  1224.  
  1225.     s[num++] = 0;
  1226.     } else {
  1227.     /* solve the resolvent cubic ... */
  1228.  
  1229.     coeffs[0] = 1.0 / 2 * r * p - 1.0 / 8 * q * q;
  1230.     coeffs[1] = -r;
  1231.     coeffs[2] = -1.0 / 2 * p;
  1232.     coeffs[3] = 1;
  1233.  
  1234.     (void) solve_cubic(coeffs, s);
  1235.  
  1236.     /* ... and take the one real solution ... */
  1237.  
  1238.     z = s[0];
  1239.  
  1240.     /* ... to build two quadric equations */
  1241.  
  1242.     u = z * z - r;
  1243.     v = 2 * z - p;
  1244.  
  1245.     if (ISZERO(u))
  1246.         u = 0;
  1247.     else if (u > 0)
  1248.         u = sqrt(u);
  1249.     else
  1250.         return 0;
  1251.  
  1252.     if (ISZERO(v))
  1253.         v = 0;
  1254.     else if (v > 0)
  1255.         v = sqrt(v);
  1256.     else
  1257.         return 0;
  1258.  
  1259.     coeffs[0] = z - u;
  1260.     coeffs[1] = q < 0 ? -v : v;
  1261.     coeffs[2] = 1;
  1262.  
  1263.     num = solve_quadric(coeffs, s);
  1264.  
  1265.     coeffs[0] = z + u;
  1266.     coeffs[1] = q < 0 ? v : -v;
  1267.     coeffs[2] = 1;
  1268.  
  1269.     num += solve_quadric(coeffs, s + num);
  1270.     }
  1271.  
  1272.     /* resubstitute */
  1273.  
  1274.     sub = 1.0 / 4 * A;
  1275.  
  1276.     for (i = 0; i < num; ++i)
  1277.     s[i] -= sub;
  1278.  
  1279.     return num;
  1280. }
  1281.  
  1282.  
  1283. int
  1284. dcomp(const void *p, const void *q)
  1285. {
  1286.     double          a,
  1287.                     b;
  1288.  
  1289.     a = *(double *) p;
  1290.     b = *(double *) q;
  1291.     if (a < b)
  1292.     return -1;
  1293.     if (a > b)
  1294.     return +1;
  1295.     return 0;
  1296. }
  1297.  
  1298. PRIVATE int
  1299. solve_rt_quartic(double *coef, double *roots)
  1300. {
  1301.     int             n,
  1302.                     i,
  1303.                     j;
  1304.     double          r[4];
  1305.  
  1306.     n = solve_quartic(coef, r);
  1307.     qsort(r, n, sizeof(double), dcomp);
  1308.  
  1309.     j = 0;
  1310.     for (i = 0; i < n; i++)
  1311.     if (r[i] > tolerance)
  1312.         roots[j++] = r[i];
  1313.     return j;
  1314. }
  1315.  
  1316. PRIVATE int
  1317. solve_rt_cubic(double *coeffs, double *roots)
  1318. {
  1319.     int             i;
  1320.  
  1321.     int             n,
  1322.                     j;
  1323.     double          r[4];
  1324.  
  1325.     n = solve_cubic(coeffs, r);
  1326.     qsort(r, n, sizeof(double), dcomp);
  1327.  
  1328.     j = 0;
  1329.     for (i = 0; i < n; i++)
  1330.     if (r[i] > tolerance)
  1331.         roots[j++] = r[i];
  1332.     return j;
  1333. }
  1334.  
  1335. /*
  1336.  * the debugging routines are still here, if you need them.
  1337.  */
  1338. #ifdef DEBUG
  1339.  
  1340. /* print solutions, and check wether they are solutions */
  1341. PUBLIC void
  1342. print_sols(double *s, int nosols, spoly p)
  1343. {
  1344.  
  1345.     int             i;
  1346.  
  1347.     if (nosols <= 0) {
  1348.     printf("no solutions\n");
  1349.     return;
  1350.     }
  1351.     printf("%d solution(s)\n", nosols);
  1352.  
  1353.     for (i = 0; i < nosols; i++)
  1354.     printf("%d: p(%lf) = %lf\n", i + 1, s[i], poly_eval(s[i], p.ord, p.coef));
  1355.  
  1356.     printf("\n");
  1357. }
  1358.  
  1359. #endif
  1360.