home *** CD-ROM | disk | FTP | other *** search
/ Otherware / Otherware_1_SB_Development.iso / mac / developm / source / macraysh.sit / Code / Source / roots.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-11-08  |  4.6 KB  |  247 lines
  1. /*
  2.  *  Roots3And4.c
  3.  *
  4.  *  Utility functions to find cubic and quartic roots,
  5.  *  coefficients are passed like this:
  6.  *
  7.  *      c[0] + c[1]*x + c[2]*x^2 + c[3]*x^3 + c[4]*x^4 = 0
  8.  *
  9.  *  The functions return the number of non-complex roots and
  10.  *  put the values into the s array.
  11.  *
  12.  *  Author:         Jochen Schwarze (schwarze@isa.de)
  13.  *
  14.  *  Jan 26, 1990    Version for Graphics Gems
  15.  *  Oct 11, 1990    Fixed sign problem for negative q's in SolveQuartic
  16.  *                  (reported by Mark Podlipec),
  17.  *                  Old-style function definitions,
  18.  *                  IsZero() as a macro
  19.  *  Nov 23, 1990    Some systems do not declare acos() and cbrt() in
  20.  *                  <math.h>, though the functions exist in the library.
  21.  *                  If large coefficients are used, EQN_EPS should be
  22.  *                  reduced considerably (e.g. to 1E-30), results will be
  23.  *                  correct but multiple roots might be reported more
  24.  *                  than once.
  25.  */
  26.  
  27. #include "common.h"
  28.  
  29. /* epsilon surrounding for near zero values */
  30.  
  31. /*
  32.  * In case M_PI isn't defined in math.h
  33.  */
  34. #ifndef M_PI
  35. #define M_PI PI
  36. #endif
  37.  
  38. #define     EQN_EPS     1e-9
  39. #define        IsZero(x)    ((x) > -EQN_EPS && (x) < EQN_EPS)
  40.  
  41. #ifndef CBRT
  42. #define     cbrt(x)     ((x) > 0.0 ? pow((double)(x), 1.0/3.0) : \
  43.               ((x) < 0.0 ? -pow((double)-(x), 1.0/3.0) : 0.0))
  44. #endif
  45.  
  46.  
  47. int SolveQuadric(c, s)
  48.     double c[ 3 ];
  49.     double s[ 2 ];
  50. {
  51.     double p, q, D;
  52.  
  53.     /* normal form: x^2 + px + q = 0 */
  54.  
  55.     p = c[ 1 ] / (2 * c[ 2 ]);
  56.     q = c[ 0 ] / c[ 2 ];
  57.  
  58.     D = p * p - q;
  59.  
  60.     if (IsZero(D))
  61.     {
  62.     s[ 0 ] = - p;
  63.     return 1;
  64.     }
  65.     else if (D > 0)
  66.     {
  67.     double sqrt_D = sqrt(D);
  68.  
  69.     s[ 0 ] =   sqrt_D - p;
  70.     s[ 1 ] = - sqrt_D - p;
  71.     return 2;
  72.     }
  73.     else /* if (D < 0) */
  74.         return 0;
  75. }
  76.  
  77.  
  78. int SolveCubic(c, s)
  79.     double c[ 4 ];
  80.     double s[ 3 ];
  81. {
  82.     int     i, num;
  83.     double  sub;
  84.     double  A, B, C;
  85.     double  sq_A, p, q;
  86.     double  cb_p, D;
  87.  
  88.     /* normal form: x^3 + Ax^2 + Bx + C = 0 */
  89.  
  90.     A = c[ 2 ] / c[ 3 ];
  91.     B = c[ 1 ] / c[ 3 ];
  92.     C = c[ 0 ] / c[ 3 ];
  93.  
  94.     /*  substitute x = y - A/3 to eliminate quadric term:
  95.     x^3 +px + q = 0 */
  96.  
  97.     sq_A = A * A;
  98.     p = 1.0/3 * (- 1.0/3 * sq_A + B);
  99.     q = 1.0/2 * (2.0/27 * A * sq_A - 1.0/3 * A * B + C);
  100.  
  101.     /* use Cardano's formula */
  102.  
  103.     cb_p = p * p * p;
  104.     D = q * q + cb_p;
  105.  
  106.     if (IsZero(D))
  107.     {
  108.     if (IsZero(q)) /* one triple solution */
  109.     {
  110.         s[ 0 ] = 0;
  111.         num = 1;
  112.     }
  113.     else /* one single and one double solution */
  114.     {
  115.         double u = cbrt(-q);
  116.         s[ 0 ] = 2 * u;
  117.         s[ 1 ] = - u;
  118.         num = 2;
  119.     }
  120.     }
  121.     else if (D < 0) /* Casus irreducibilis: three real solutions */
  122.     {
  123.     double phi = 1.0/3 * acos(-q / sqrt(-cb_p));
  124.     double t = 2 * sqrt(-p);
  125.  
  126.     s[ 0 ] =   t * cos(phi);
  127.     s[ 1 ] = - t * cos(phi + M_PI / 3);
  128.     s[ 2 ] = - t * cos(phi - M_PI / 3);
  129.     num = 3;
  130.     }
  131.     else /* one real solution */
  132.     {
  133.     double sqrt_D = sqrt(D);
  134.     double u = cbrt(sqrt_D - q);
  135.     double v = - cbrt(sqrt_D + q);
  136.  
  137.     s[ 0 ] = u + v;
  138.     num = 1;
  139.     }
  140.  
  141.     /* resubstitute */
  142.  
  143.     sub = 1.0/3 * A;
  144.  
  145.     for (i = 0; i < num; ++i)
  146.     s[ i ] -= sub;
  147.  
  148.     return num;
  149. }
  150.  
  151.  
  152. int SolveQuartic(c, s)
  153.     double c[ 5 ]; 
  154.     double s[ 4 ];
  155. {
  156.     double  coeffs[ 4 ];
  157.     double  z, u, v, sub;
  158.     double  A, B, C, D;
  159.     double  sq_A, p, q, r;
  160.     int     i, num;
  161.  
  162.     /* normal form: x^4 + Ax^3 + Bx^2 + Cx + D = 0 */
  163.  
  164.     A = c[ 3 ] / c[ 4 ];
  165.     B = c[ 2 ] / c[ 4 ];
  166.     C = c[ 1 ] / c[ 4 ];
  167.     D = c[ 0 ] / c[ 4 ];
  168.  
  169.     /*  substitute x = y - A/4 to eliminate cubic term:
  170.     x^4 + px^2 + qx + r = 0 */
  171.  
  172.     sq_A = A * A;
  173.     p = - 3.0/8 * sq_A + B;
  174.     q = 1.0/8 * sq_A * A - 1.0/2 * A * B + C;
  175.     r = - 3.0/256*sq_A*sq_A + 1.0/16*sq_A*B - 1.0/4*A*C + D;
  176.  
  177.     if (IsZero(r))
  178.     {
  179.     /* no absolute term: y(y^3 + py + q) = 0 */
  180.  
  181.     coeffs[ 0 ] = q;
  182.     coeffs[ 1 ] = p;
  183.     coeffs[ 2 ] = 0;
  184.     coeffs[ 3 ] = 1;
  185.  
  186.     num = SolveCubic(coeffs, s);
  187.  
  188.     s[ num++ ] = 0;
  189.     }
  190.     else
  191.     {
  192.     /* solve the resolvent cubic ... */
  193.  
  194.     coeffs[ 0 ] = 1.0/2 * r * p - 1.0/8 * q * q;
  195.     coeffs[ 1 ] = - r;
  196.     coeffs[ 2 ] = - 1.0/2 * p;
  197.     coeffs[ 3 ] = 1;
  198.  
  199.     (void) SolveCubic(coeffs, s);
  200.  
  201.     /* ... and take the one real solution ... */
  202.  
  203.     z = s[ 0 ];
  204.  
  205.     /* ... to build two quadric equations */
  206.  
  207.     u = z * z - r;
  208.     v = 2 * z - p;
  209.  
  210.     if (IsZero(u))
  211.         u = 0;
  212.     else if (u > 0)
  213.         u = sqrt(u);
  214.     else
  215.         return 0;
  216.  
  217.     if (IsZero(v))
  218.         v = 0;
  219.     else if (v > 0)
  220.         v = sqrt(v);
  221.     else
  222.         return 0;
  223.  
  224.     coeffs[ 0 ] = z - u;
  225.     coeffs[ 1 ] = q < 0 ? -v : v;
  226.     coeffs[ 2 ] = 1;
  227.  
  228.     num = SolveQuadric(coeffs, s);
  229.  
  230.     coeffs[ 0 ]= z + u;
  231.     coeffs[ 1 ] = q < 0 ? v : -v;
  232.     coeffs[ 2 ] = 1;
  233.  
  234.     num += SolveQuadric(coeffs, s + num);
  235.     }
  236.  
  237.     /* resubstitute */
  238.  
  239.     sub = 1.0/4 * A;
  240.  
  241.     for (i = 0; i < num; ++i)
  242.     s[ i ] -= sub;
  243.  
  244.     return num;
  245. }
  246.  
  247.