| /* +++Date last modified: 05-Jul-1997 */ |
| |
| /* |
| ** CUBIC.C - Solve a cubic polynomial |
| ** public domain by Ross Cottrell |
| */ |
| |
| #include <math.h> |
| #include <stdlib.h> |
| #include "snipmath.h" |
| |
| void SolveCubic(double a, |
| double b, |
| double c, |
| double d, |
| int *solutions, |
| double *x) |
| { |
| long double a1 = b/a, a2 = c/a, a3 = d/a; |
| long double Q = (a1*a1 - 3.0*a2)/9.0; |
| long double R = (2.0*a1*a1*a1 - 9.0*a1*a2 + 27.0*a3)/54.0; |
| double R2_Q3 = R*R - Q*Q*Q; |
| |
| double theta; |
| |
| if (R2_Q3 <= 0) |
| { |
| *solutions = 3; |
| theta = acos(R/sqrt(Q*Q*Q)); |
| x[0] = -2.0*sqrt(Q)*cos(theta/3.0) - a1/3.0; |
| x[1] = -2.0*sqrt(Q)*cos((theta+2.0*PI)/3.0) - a1/3.0; |
| x[2] = -2.0*sqrt(Q)*cos((theta+4.0*PI)/3.0) - a1/3.0; |
| } |
| else |
| { |
| *solutions = 1; |
| x[0] = pow(sqrt(R2_Q3)+fabs(R), 1/3.0); |
| x[0] += Q/x[0]; |
| x[0] *= (R < 0.0) ? 1 : -1; |
| x[0] -= a1/3.0; |
| } |
| } |
| |
| #ifdef TEST |
| |
| int main(void) |
| { |
| double a1 = 1.0, b1 = -10.5, c1 = 32.0, d1 = -30.0; |
| double a2 = 1.0, b2 = -4.5, c2 = 17.0, d2 = -30.0; |
| double x[3]; |
| int solutions; |
| |
| SolveCubic(a1, b1, c1, d1, &solutions, x); |
| |
| /* should get 3 solutions: 2, 6 & 2.5 */ |
| |
| SolveCubic(a2, b2, c2, d2, &solutions, x); |
| |
| /* should get 1 solution: 2.5 */ |
| |
| return 0; |
| } |
| |
| #endif /* TEST */ |
