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