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