| 32 | } |
| 33 | |
| 34 | static int solveCubicNormed(double x[3], double a, double b, double c) { |
| 35 | double a2 = a*a; |
| 36 | double q = 1/9.*(a2-3*b); |
| 37 | double r = 1/54.*(a*(2*a2-9*b)+27*c); |
| 38 | double r2 = r*r; |
| 39 | double q3 = q*q*q; |
| 40 | a *= 1/3.; |
| 41 | if (r2 < q3) { |
| 42 | double t = r/sqrt(q3); |
| 43 | if (t < -1) t = -1; |
| 44 | if (t > 1) t = 1; |
| 45 | t = acos(t); |
| 46 | q = -2*sqrt(q); |
| 47 | x[0] = q*cos(1/3.*t)-a; |
| 48 | x[1] = q*cos(1/3.*(t+2*M_PI))-a; |
| 49 | x[2] = q*cos(1/3.*(t-2*M_PI))-a; |
| 50 | return 3; |
| 51 | } else { |
| 52 | double u = (r < 0 ? 1 : -1)*pow(fabs(r)+sqrt(r2-q3), 1/3.); |
| 53 | double v = u == 0 ? 0 : q/u; |
| 54 | x[0] = (u+v)-a; |
| 55 | if (u == v || fabs(u-v) < 1e-12*fabs(u+v)) { |
| 56 | x[1] = -.5*(u+v)-a; |
| 57 | return 2; |
| 58 | } |
| 59 | return 1; |
| 60 | } |
| 61 | } |
| 62 | |
| 63 | int solveCubic(double x[3], double a, double b, double c, double d) { |
| 64 | if (a != 0) { |