| |
| #include <cmath> |
| #include <iostream> |
| |
| #include "real.h" |
| |
| |
| int solve_equation(Real* roots, const Real& k0, const Real& k1) { |
| if (equal(k1, 0)) return 0; |
| if (roots) roots[0] = -k0/k1; |
| return 1; |
| } |
| |
| int solve_equation(Real* roots, const Real& k0, const Real& k1, const Real& k2) { |
| if (equal(k2, 0)) return solve_equation(roots,k0, k1); |
| Real D = k1*k1 - 4*k2*k0; |
| if (equal(D, 0)) { |
| if (roots) roots[0] = -0.5*k1/k2; |
| return 1; |
| } else |
| if (D > 0) { |
| if (roots) { |
| Real a = sqrt(D); |
| Real b = -0.5/k2; |
| roots[0] = (k1 - a)*b; |
| roots[1] = (k1 + a)*b; |
| } |
| return 2; |
| } |
| return 0; |
| } |
| |
| int solve_equation(Real* roots, const Real& k0, const Real& k1, const Real& k2, const Real& k3) { |
| if (equal(k3, 0)) return solve_equation(roots, k0, k1, k2); |
| |
| Real k = 1/k3; |
| Real a = k2*k; |
| Real b = k1*k; |
| Real c = k0*k; |
| |
| Real Q = (a*a - 3*b)/9; |
| Real Q3 = Q*Q*Q; |
| Real R = (2*a*a*a - 9*a*b + 27*c)/54; |
| Real S = Q3 - R*R; |
| |
| if (equal(S, 0)) { |
| if (roots) { |
| Real rr = cbrt(R); |
| Real aa = -a/3; |
| roots[0] = aa - 2*rr; |
| roots[1] = aa + rr; |
| } |
| return 2; |
| } else |
| if (S > 0) { |
| if (roots) { |
| Real ph = acos(R/sqrt(Q3))/3; |
| Real qq = -2*sqrt(Q); |
| Real aa = -a/3; |
| roots[0] = qq*cos(ph) + aa; |
| roots[1] = qq*cos(ph + 2*pi/3) + aa; |
| roots[2] = qq*cos(ph - 2*pi/3) + aa; |
| } |
| return 3; |
| } else |
| if (equal(Q, 0)) { |
| if (roots) roots[0] = -cbrt(c - a*a*a/27) - a/3; |
| return 1; |
| } else |
| if (Q > 0) { |
| if (roots) { |
| Real ph = acosh( fabs(R)/sqrt(Q3) )/3; |
| roots[0] = -2*sign(R)*sqrt(Q)*cosh(ph) - a/3; |
| } |
| return 1; |
| } else { |
| if (roots) { |
| Real ph = asinh( fabs(R)/sqrt(-Q3) )/3; |
| roots[0] = -2*sign(R)*sqrt(-Q)*sinh(ph) - a/3; |
| } |
| return 1; |
| } |
| return 0; |
| } |