home *** CD-ROM | disk | FTP | other *** search
/ Amiga Tools 3 / Amiga Tools 3.iso / grafik / raytracing / rayshade-4.0.6.3 / libray / libobj / roots.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-09  |  4.8 KB  |  254 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 "libcommon/common.h"
  28.  
  29. #ifdef AMIGA
  30. extern double    cbrt();
  31. #else
  32. extern double   sqrt(), cbrt(), cos(), acos();
  33. #endif
  34.  
  35. /* epsilon surrounding for near zero values */
  36.  
  37. /*
  38.  * In case M_PI isn't defined in math.h
  39.  */
  40. #ifndef M_PI
  41. #define M_PI PI
  42. #endif
  43.  
  44. #define     EQN_EPS     1e-9
  45. #define        IsZero(x)    ((x) > -EQN_EPS && (x) < EQN_EPS)
  46. #define        IsZero2(x)    ((x) > -EQN_EPS*EQN_EPS && (x) < EQN_EPS*EQN_EPS)
  47.  
  48. #ifndef CBRT
  49. #define     cbrt(x)     ((x) > 0.0 ? pow((double)(x), 1.0/3.0) : \
  50.               ((x) < 0.0 ? -pow((double)-(x), 1.0/3.0) : 0.0))
  51. #endif
  52.  
  53.  
  54. int SolveQuadric(c, s)
  55.     double c[ 3 ];
  56.     double s[ 2 ];
  57. {
  58.     double p, q, D;
  59.  
  60.     /* normal form: x^2 + px + q = 0 */
  61.  
  62.     p = c[ 1 ] / (2 * c[ 2 ]);
  63.     q = c[ 0 ] / c[ 2 ];
  64.  
  65.     D = p * p - q;
  66.  
  67.     if (IsZero(D))
  68.     {
  69.     s[ 0 ] = - p;
  70.     return 1;
  71.     }
  72.     else if (D > 0)
  73.     {
  74.     double sqrt_D = sqrt(D);
  75.  
  76.     s[ 0 ] =   sqrt_D - p;
  77.     s[ 1 ] = - sqrt_D - p;
  78.     return 2;
  79.     }
  80.     else /* if (D < 0) */
  81.         return 0;
  82. }
  83.  
  84.  
  85. int SolveCubic(c, s)
  86.     double c[ 4 ];
  87.     double s[ 3 ];
  88. {
  89.     int     i, num;
  90.     double  sub;
  91.     double  A, B, C;
  92.     double  sq_A, p, q;
  93.     double  cb_p, D;
  94.  
  95.     /* normal form: x^3 + Ax^2 + Bx + C = 0 */
  96.  
  97.     A = c[ 2 ] / c[ 3 ];
  98.     B = c[ 1 ] / c[ 3 ];
  99.     C = c[ 0 ] / c[ 3 ];
  100.  
  101.     /*  substitute x = y - A/3 to eliminate quadric term:
  102.     x^3 +px + q = 0 */
  103.  
  104.     sq_A = A * A;
  105.     p = 1.0/3 * (- 1.0/3 * sq_A + B);
  106.     q = 1.0/2 * (2.0/27 * A * sq_A - 1.0/3 * A * B + C);
  107.  
  108.     /* use Cardano's formula */
  109.  
  110.     cb_p = p * p * p;
  111.     D = q * q + cb_p;
  112.  
  113.     if (IsZero2(D))
  114.     {
  115.     if (IsZero(q)) /* one triple solution */
  116.     {
  117.         s[ 0 ] = 0;
  118.         num = 1;
  119.     }
  120.     else /* one single and one double solution */
  121.     {
  122.         double u = cbrt(-q);
  123.         s[ 0 ] = 2 * u;
  124.         s[ 1 ] = - u;
  125.         num = 2;
  126.     }
  127.     }
  128.     else if (D < 0) /* Casus irreducibilis: three real solutions */
  129.     {
  130.     double phi = 1.0/3 * acos(-q / sqrt(-cb_p));
  131.     double t = 2 * sqrt(-p);
  132.  
  133.     s[ 0 ] =   t * cos(phi);
  134.     s[ 1 ] = - t * cos(phi + M_PI / 3);
  135.     s[ 2 ] = - t * cos(phi - M_PI / 3);
  136.     num = 3;
  137.     }
  138.     else /* one real solution */
  139.     {
  140.     double sqrt_D = sqrt(D);
  141.     double u = cbrt(sqrt_D - q);
  142.     double v = - cbrt(sqrt_D + q);
  143.  
  144.     s[ 0 ] = u + v;
  145.     num = 1;
  146.     }
  147.  
  148.     /* resubstitute */
  149.  
  150.     sub = 1.0/3 * A;
  151.  
  152.     for (i = 0; i < num; ++i)
  153.     s[ i ] -= sub;
  154.  
  155.     return num;
  156. }
  157.  
  158.  
  159. int SolveQuartic(c, s)
  160.     double c[ 5 ]; 
  161.     double s[ 4 ];
  162. {
  163.     double  coeffs[ 4 ];
  164.     double  z, u, v, sub;
  165.     double  A, B, C, D;
  166.     double  sq_A, p, q, r;
  167.     int     i, num;
  168.  
  169.     /* normal form: x^4 + Ax^3 + Bx^2 + Cx + D = 0 */
  170.  
  171.     A = c[ 3 ] / c[ 4 ];
  172.     B = c[ 2 ] / c[ 4 ];
  173.     C = c[ 1 ] / c[ 4 ];
  174.     D = c[ 0 ] / c[ 4 ];
  175.  
  176.     /*  substitute x = y - A/4 to eliminate cubic term:
  177.     x^4 + px^2 + qx + r = 0 */
  178.  
  179.     sq_A = A * A;
  180.     p = - 3.0/8 * sq_A + B;
  181.     q = 1.0/8 * sq_A * A - 1.0/2 * A * B + C;
  182.     r = - 3.0/256*sq_A*sq_A + 1.0/16*sq_A*B - 1.0/4*A*C + D;
  183.  
  184.     if (IsZero(r))
  185.     {
  186.     /* no absolute term: y(y^3 + py + q) = 0 */
  187.  
  188.     coeffs[ 0 ] = q;
  189.     coeffs[ 1 ] = p;
  190.     coeffs[ 2 ] = 0;
  191.     coeffs[ 3 ] = 1;
  192.  
  193.     num = SolveCubic(coeffs, s);
  194.  
  195.     s[ num++ ] = 0;
  196.     }
  197.     else
  198.     {
  199.     /* solve the resolvent cubic ... */
  200.  
  201.     coeffs[ 0 ] = 1.0/2 * r * p - 1.0/8 * q * q;
  202.     coeffs[ 1 ] = - r;
  203.     coeffs[ 2 ] = - 1.0/2 * p;
  204.     coeffs[ 3 ] = 1;
  205.  
  206.     (void) SolveCubic(coeffs, s);
  207.  
  208.     /* ... and take the one real solution ... */
  209.  
  210.     z = s[ 0 ];
  211.  
  212.     /* ... to build two quadric equations */
  213.  
  214.     u = z * z - r;
  215.     v = 2 * z - p;
  216.  
  217.     if (IsZero(u))
  218.         u = 0;
  219.     else if (u > 0)
  220.         u = sqrt(u);
  221.     else
  222.         return 0;
  223.  
  224.     if (IsZero(v))
  225.         v = 0;
  226.     else if (v > 0)
  227.         v = sqrt(v);
  228.     else
  229.         return 0;
  230.  
  231.     coeffs[ 0 ] = z - u;
  232.     coeffs[ 1 ] = q < 0 ? -v : v;
  233.     coeffs[ 2 ] = 1;
  234.  
  235.     num = SolveQuadric(coeffs, s);
  236.  
  237.     coeffs[ 0 ]= z + u;
  238.     coeffs[ 1 ] = q < 0 ? v : -v;
  239.     coeffs[ 2 ] = 1;
  240.  
  241.     num += SolveQuadric(coeffs, s + num);
  242.     }
  243.  
  244.     /* resubstitute */
  245.  
  246.     sub = 1.0/4 * A;
  247.  
  248.     for (i = 0; i < num; ++i)
  249.     s[ i ] -= sub;
  250.  
  251.     return num;
  252. }
  253.  
  254.