27 #ifndef OPM_POLYNOMIAL_UTILS_HH
28 #define OPM_POLYNOMIAL_UTILS_HH
50 template <
class Scalar,
class SolContainer>
55 if (std::abs(scalarValue(a)) < 1e-30)
78 template <
class Scalar,
class SolContainer>
85 if (std::abs(scalarValue(a)) < 1e-30)
89 Scalar Delta = b*b - 4*a*c;
94 sol[0] = (- b + Delta)/(2*a);
95 sol[1] = (- b - Delta)/(2*a);
99 std::swap(sol[0], sol[1]);
104 template <
class Scalar,
class SolContainer>
105 void invertCubicPolynomialPostProcess_(SolContainer& sol,
114 for (
int i = 0; i < numSol; ++i) {
116 Scalar fOld = d + x*(c + x*(b + x*a));
118 Scalar fPrime = c + x*(2*b + x*3*a);
119 if (std::abs(scalarValue(fPrime)) < 1e-30)
123 Scalar fNew = d + x*(c + x*(b + x*a));
124 if (std::abs(scalarValue(fNew)) < std::abs(scalarValue(fOld)))
147 template <
class Scalar,
class SolContainer>
155 if (std::abs(scalarValue(a)) < 1e-30)
165 Scalar p = c - b*b/3;
166 Scalar q = d + (2*b*b*b - 9*b*c)/27;
168 if (std::abs(scalarValue(p)) > 1e-30 && std::abs(scalarValue(q)) > 1e-30) {
199 Scalar wDisc = q*q/4 + p*p*p/27;
202 Scalar u = - q/2 + sqrt(wDisc);
203 if (u < 0) u = - pow(-u, 1.0/3);
204 else u = pow(u, 1.0/3);
208 sol[0] = u - p/(3*u) - b/3;
212 invertCubicPolynomialPostProcess_(sol, 1, a, b, c, d);
216 Scalar uCubedRe = - q/2;
217 Scalar uCubedIm = sqrt(-wDisc);
219 Scalar uAbs = pow(sqrt(uCubedRe*uCubedRe + uCubedIm*uCubedIm), 1.0/3);
220 Scalar phi = atan2(uCubedIm, uCubedRe)/3;
261 for (
int i = 0; i < 3; ++i) {
262 sol[i] = cos(phi)*(uAbs - p/(3*uAbs)) - b/3;
268 invertCubicPolynomialPostProcess_(sol, 3, a, b, c, d);
271 std::sort(sol, sol + 3);
278 else if (std::abs(scalarValue(p)) < 1e-30 && std::abs(scalarValue(q)) < 1e-30) {
280 sol[0] = sol[1] = sol[2] = 0.0 - b/3;
283 else if (std::abs(scalarValue(p)) < 1e-30 && std::abs(scalarValue(q)) > 1e-30) {
288 if (-q > 0) t = pow(-q, 1./3);
289 else t = - pow(q, 1./3);
295 assert(std::abs(scalarValue(p)) > 1e-30 && std::abs(scalarValue(q)) > 1e-30);
306 sol[0] = -sqrt(-p) - b/3;;
308 sol[2] = sqrt(-p) - b/3;
330 template <
class Scalar,
class SolContainer>
338 if (std::abs(scalarValue(a)) < 1e-30)
345 Scalar p = (3.0 * a * c - b * b) / (3.0 * a * a);
346 Scalar q = (2.0 * b * b * b - 9.0 * a * b * c + 27.0 * d * a * a) / (27.0 * a * a * a);
350 Scalar discr = 4.0 * p * p * p + 27.0 * q * q;
354 Scalar theta = (1.0 / 3.0) * acos( ((3.0 * q) / (2.0 * p)) * sqrt(-3.0 / p) );
357 sol[0] = 2.0 * sqrt(-p / 3.0) * cos( theta ) - b / (3.0 * a);
358 sol[1] = 2.0 * sqrt(-p / 3.0) * cos( theta - ((2.0 * M_PI) / 3.0) ) - b / (3.0 * a);
359 sol[2] = 2.0 * sqrt(-p / 3.0) * cos( theta - ((4.0 * M_PI) / 3.0) ) - b / (3.0 * a);
362 std::sort(sol, sol + 3);
368 else if (discr > 0.0) {
375 Scalar theta = (1.0 / 3.0) * acosh( ((-3.0 * abs(q)) / (2.0 * p)) * sqrt(-3.0 / p) );
378 t = ( (-2.0 * abs(q)) / q ) * sqrt(-p / 3.0) * cosh(theta);
382 Scalar theta = (1.0 / 3.0) * asinh( ((3.0 * q) / (2.0 * p)) * sqrt(3.0 / p) );
384 t = -2.0 * sqrt(p / 3.0) * sinh(theta);
388 std::runtime_error(
" p = 0 in cubic root solver!");
392 sol[0] = t - b / (3.0 * a);
400 sol[0] = sol[1] = sol[2] = 0.0 - b / (3.0 * a);
404 sol[0] = (3.0 * q / p) - b / (3.0 * a);
405 sol[1] = sol[2] = (-3.0 * q) / (2.0 * p) - b / (3.0 * a);
406 std::sort(sol, sol + 3);
unsigned invertQuadraticPolynomial(SolContainer &sol, Scalar a, Scalar b, Scalar c)
Invert a quadratic polynomial analytically.
Definition: PolynomialUtils.hpp:79
unsigned invertLinearPolynomial(SolContainer &sol, Scalar a, Scalar b)
Invert a linear polynomial analytically.
Definition: PolynomialUtils.hpp:51
unsigned cubicRoots(SolContainer *sol, Scalar a, Scalar b, Scalar c, Scalar d)
Invert a cubic polynomial analytically.
Definition: PolynomialUtils.hpp:331
unsigned invertCubicPolynomial(SolContainer *sol, Scalar a, Scalar b, Scalar c, Scalar d)
Invert a cubic polynomial analytically.
Definition: PolynomialUtils.hpp:148