home *** CD-ROM | disk | FTP | other *** search
/ Point Programming 1 / PPROG1.ISO / c / snippets / cubic.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-04-03  |  1.5 KB  |  62 lines
  1. /*
  2. **  CUBIC.C - Solve a cubic polynomial
  3. **  public domain by Ross Cottrell
  4. */
  5.  
  6. #include <math.h>
  7. #include <stdlib.h>
  8.  
  9. void SolveCubic(double  a,
  10.                 double  b,
  11.                 double  c,
  12.                 double  d,
  13.                 int    *solutions,
  14.                 double *x)
  15. {
  16.       long double    a1 = b/a, a2 = c/a, a3 = d/a;
  17.       long double    Q = (a1*a1 - 3.0*a2)/9.0;
  18.       long double R = (2.0*a1*a1*a1 - 9.0*a1*a2 + 27.0*a3)/54.0;
  19.       double    R2_Q3 = R*R - Q*Q*Q;
  20.  
  21.       double    theta;
  22.  
  23.       if (R2_Q3 <= 0)
  24.       {
  25.             *solutions = 3;
  26.             theta = acos(R/sqrt(Q*Q*Q));
  27.             x[0] = -2.0*sqrt(Q)*cos(theta/3.0) - a1/3.0;
  28.             x[1] = -2.0*sqrt(Q)*cos((theta+2.0*M_PI)/3.0) - a1/3.0;
  29.             x[2] = -2.0*sqrt(Q)*cos((theta+4.0*M_PI)/3.0) - a1/3.0;
  30.       }
  31.       else
  32.       {
  33.             *solutions = 1;
  34.             x[0] = pow(sqrt(R2_Q3)+fabs(R), 1/3.0);
  35.             x[0] += Q/x[0];
  36.             x[0] *= (R < 0.0) ? 1 : -1;
  37.             x[0] -= a1/3.0;
  38.       }
  39. }
  40.  
  41. #ifdef TEST
  42.  
  43. int main(void)
  44. {
  45.       double  a1 = 1.0, b1 = -10.5, c1 = 32.0, d1 = -30.0;
  46.       double  a2 = 1.0, b2 = -4.5, c2 = 17.0, d2 = -30.0;
  47.       double  x[3];
  48.       int     solutions;
  49.  
  50.       SolveCubic(a1, b1, c1, d1, &solutions, x);
  51.  
  52.       /* should get 3 solutions: 2, 6 & 2.5   */
  53.  
  54.       SolveCubic(a2, b2, c2, d2, &solutions, x);
  55.  
  56.       /* should get 1 solution: 2.5           */
  57.  
  58.       return 0;
  59. }
  60.  
  61. #endif /* TEST */
  62.