home *** CD-ROM | disk | FTP | other *** search
/ Photo CD Demo 1 / Demo.bin / gems / graphics / roots3an.c < prev    next >
C/C++ Source or Header  |  1992-04-09  |  4KB  |  235 lines
  1. /*
  2. Cubic and Quartic Roots
  3. by Jochen Schwarze
  4. from "Graphics Gems", Academic Press, 1990
  5. */
  6.  
  7. /*
  8. **  Utility functions to find cubic and quartic roots,
  9. **  parameters are given like this:
  10. **
  11. **      c[0] + c[1]*x + c[2]*x^2 + c[3]*x^3 + c[4]*x^4 = 0
  12. **
  13. **  The functions return the number of non-complex roots and
  14. **  put the values into the s array.
  15. **
  16. **  Author:         Jochen Schwarze (schwarze@isaak.isa.de)
  17. **  Last modified:  Jan 26, 1990
  18. */
  19.  
  20. #include    <math.h>
  21.  
  22. /* epsilon surrounding for near zero values */
  23.  
  24. #define     EQN_EPS     1e-9
  25.  
  26. static int
  27. IsZero(double x)
  28. {
  29.     return x > - EQN_EPS && x < EQN_EPS;
  30. }
  31.  
  32.  
  33. int
  34. SolveQuadric(double c[ 3 ], double s[ 2 ])
  35. {
  36.     double p, q, D;
  37.  
  38.     /* normal form: x^2 + px + q = 0 */
  39.  
  40.     p = c[ 1 ] / (2 * c[ 2 ]);
  41.     q = c[ 0 ] / c[ 2 ];
  42.  
  43.  
  44.     D = p * p - q;
  45.  
  46.     if     (IsZero(D))
  47.     {
  48.         s[ 0 ] = - p;
  49.         return 1;
  50.     }
  51.     else if (D < 0)
  52.     {
  53.         return 0;
  54.     }
  55.     else if (D > 0)
  56.     {
  57.         double sqrt_D = sqrt(D);
  58.  
  59.         s[ 0 ] =   sqrt_D - p;
  60.         s[ 1 ] = - sqrt_D - p;
  61.         return 2;
  62.     }
  63. }
  64.  
  65.  
  66. int
  67. SolveCubic(double c[ 4 ], double s[ 3 ])
  68. {
  69.     int     i, num;
  70.     double  sub;
  71.     double  A, B, C;
  72.     double  sq_A, p, q;
  73.     double  cb_p, D;
  74.  
  75.     /* normal form: x^3 + Ax^2 + Bx + C = 0 */
  76.  
  77.     A = c[ 2 ] / c[ 3 ];
  78.     B = c[ 1 ] / c[ 3 ];
  79.     C = c[ 0 ] / c[ 3 ];
  80.  
  81.     /*  substitute x = y - A/3 to eliminate quadric term:
  82.     x^3 +px + q = 0 */
  83.  
  84.     sq_A = A * A;
  85.     p = 1.0/3 * (- 1.0/3 * sq_A + B);
  86.     q = 1.0/2 * (2.0/27 * A * sq_A - 1.0/3 * A * B + C);
  87.  
  88.  
  89.     /* use Cardano's formula */
  90.  
  91.     cb_p = p * p * p;
  92.     D = q * q + cb_p;
  93.  
  94.     if     (IsZero(D))
  95.     {
  96.         if (IsZero(q))     /* one triple solution */
  97.         {
  98.             s[ 0 ] = 0;
  99.             num = 1;
  100.         }
  101.         else         /* one single and one double solution */
  102.         {
  103.             double u = cbrt(-q);
  104.             s[ 0 ] = 2 * u;
  105.             s[ 1 ] = - u;
  106.             num = 2;
  107.         }
  108.     }
  109.     else if (D < 0) /* Casus irreducibilis: three real solutions */
  110.     {
  111.         double phi = 1.0/3 * acos(-q / sqrt(-cb_p));
  112.         double t = 2 * sqrt(-p);
  113.  
  114.         s[ 0 ] =   t * cos(phi);
  115.         s[ 1 ] = - t * cos(phi + M_PI / 3);
  116.         s[ 2 ] = - t * cos(phi - M_PI / 3);
  117.         num = 3;
  118.     }
  119.     else             /* one real solution */
  120.     {
  121.         double sqrt_D = sqrt(D);
  122.         double u = cbrt(sqrt_D - q);
  123.         double v = - cbrt(sqrt_D + q);
  124.  
  125.         s[ 0 ] = u + v;
  126.         num = 1;
  127.     }
  128.  
  129.     /* resubstitute */
  130.  
  131.     sub = 1.0/3 * A;
  132.  
  133.     for (i = 0; i < num; ++i)
  134.         s[ i ] -= sub;
  135.  
  136.     return num;
  137. }
  138.  
  139.  
  140.  
  141. int
  142. SolveQuartic(double c[ 5 ], double s[ 4 ])
  143. {
  144.     double  coeffs[ 4 ];
  145.     double  z, u, v, sub;
  146.     double  A, B, C, D;
  147.     double  sq_A, p, q, r;
  148.     int     i, num;
  149.  
  150.     /* normal form: x^4 + Ax^3 + Bx^2 + Cx + D = 0 */
  151.  
  152.     A = c[ 3 ] / c[ 4 ];
  153.     B = c[ 2 ] / c[ 4 ];
  154.     C = c[ 1 ] / c[ 4 ];
  155.     D = c[ 0 ] / c[ 4 ];
  156.  
  157.     /*  substitute x = y - A/4 to eliminate cubic term:
  158.     x^4 + px^2 + qx + r = 0 */
  159.  
  160.     sq_A = A * A;
  161.     p = - 3.0/8 * sq_A + B;
  162.     q = 1.0/8 * sq_A * A - 1.0/2 * A * B + C;
  163.     r = - 3.0/256*sq_A*sq_A + 1.0/16*sq_A*B - 1.0/4*A*C + D;
  164.  
  165.     if (IsZero(r))
  166.     {
  167.         /* no absolute term: y(y^3 + py + q) = 0 */
  168.  
  169.         coeffs[ 0 ] = q;
  170.         coeffs[ 1 ] = p;
  171.         coeffs[ 2 ] = 0;
  172.         coeffs[ 3 ] = 1;
  173.  
  174.         num = SolveCubic(coeffs, s);
  175.  
  176.         s[ num++ ] = 0;
  177.     }
  178.     else
  179.     {
  180.         /* solve the resolvent cubic ... */
  181.  
  182.         coeffs[ 0 ] = 1.0/2 * r * p - 1.0/8 * q * q;
  183.         coeffs[ 1 ] = - r;
  184.         coeffs[ 2 ] = - 1.0/2 * p;
  185.         coeffs[ 3 ] = 1;
  186.  
  187.         (void) SolveCubic(coeffs, s);
  188.  
  189.  
  190.         /* ... and take the one real solution ... */
  191.  
  192.         z = s[ 0 ];
  193.  
  194.         /* ... to build two quadric equations */
  195.  
  196.         u = z * z - r;
  197.         v = 2 * z - p;
  198.  
  199.         if     (IsZero(u))
  200.             u = 0;
  201.         else if (u > 0)
  202.             u = sqrt(u);
  203.         else
  204.             return 0;
  205.  
  206.         if (IsZero(v))
  207.             v = 0;
  208.         else if (v > 0)
  209.             v = sqrt(v);
  210.         else
  211.             return 0;
  212.  
  213.         coeffs[ 0 ] = z - u;
  214.         coeffs[ 1 ] = v;
  215.         coeffs[ 2 ] = 1;
  216.  
  217.         num = SolveQuadric(coeffs, s);
  218.  
  219.         coeffs[ 0 ]= z + u;
  220.         coeffs[ 1 ] = - v;
  221.         coeffs[ 2 ] = 1;
  222.  
  223.         num += SolveQuadric(coeffs, s + num);
  224.     }
  225.  
  226.         /* resubstitute */
  227.  
  228.     sub = 1.0/4 * A;
  229.  
  230.     for (i = 0; i < num; ++i)
  231.         s[ i ] -= sub;
  232.  
  233.     return num;
  234. }
  235.