/* ** CUBIC.C - Solve a cubic polynomial */ #include <math.h> #include <stdlib.h> void SolveCubic(double a, /* coefficient of x^3 */ double b, /* coefficient of x^2 */ double c, /* coefficient of x */ double d, /* constant term */ int *solutions, /* # of distinct solutions */ double *x) /* array of solutions */ { 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; } }