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