blob: 7740eb2d261bf70f9c132a3c0c81716126d45360 [file] [log] [blame]
 /* +++Date last modified: 05-Jul-1997 */ /* ** CUBIC.C - Solve a cubic polynomial ** public domain by Ross Cottrell */ #include #include #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 */