25int solve3(
double *coeff,
double *roots)
29 double p, q, disc, b_over_3a, c_over_a, d_over_a;
30 double r, theta, temp,
alpha, beta;
32 a = coeff[3], b = coeff[2], c = coeff[1], d = coeff[0];
34 return solve2(coeff, roots);
35 b_over_3a = b / (3 * a);
39 p = b_over_3a * b_over_3a;
40 q = 2 * b_over_3a * p - b_over_3a * c_over_a + d_over_a;
42 disc = q * q + 4 * p * p * p;
45 r = .5 * sqrt(-disc + q * q);
46 theta = atan2(sqrt(-disc), -q);
48 roots[0] = temp * cos(theta / 3);
49 roots[1] = temp * cos((theta +
M_PI +
M_PI) / 3);
50 roots[2] = temp * cos((theta -
M_PI -
M_PI) / 3);
53 alpha = .5 * (sqrt(disc) - q);
55 roots[0] = cbrt(
alpha) + cbrt(beta);
59 roots[1] = roots[2] = -.5 * roots[0], rootn = 3;
62 for (i = 0; i < rootn; i++)
63 roots[i] -= b_over_3a;
68static int solve2(
double *coeff,
double *roots)
71 double disc, b_over_2a, c_over_a;
73 a = coeff[2], b = coeff[1], c = coeff[0];
75 return solve1(coeff, roots);
76 b_over_2a = b / (2 * a);
79 disc = b_over_2a * b_over_2a - c_over_a;
82 }
else if (disc > 0) {
83 roots[0] = -b_over_2a + sqrt(disc);
84 roots[1] = -2 * b_over_2a - roots[0];
87 roots[0] = -b_over_2a;