27#ifndef OPM_POLYNOMIAL_UTILS_HH
28#define OPM_POLYNOMIAL_UTILS_HH
50template <
class Scalar,
class SolContainer>
78template <
class Scalar,
class SolContainer>
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]);
104template <
class Scalar,
class SolContainer>
105void 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);
123 Scalar fNew = d + x*(c + x*(b + x*a));
147template <
class Scalar,
class SolContainer>
165 Scalar p = c - b*b/3;
166 Scalar q = d + (2*b*b*b - 9*b*c)/27;
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);
280 sol[0] = sol[1] = sol[2] = 0.0 - b/3;
288 if (-q > 0) t =
pow(-q, 1./3);
289 else t = -
pow(q, 1./3);
306 sol[0] = -
sqrt(-p) - b/3;;
308 sol[2] =
sqrt(-p) - b/3;
330template <
class Scalar,
class SolContainer>
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);
Definition: Air_Mesitylene.hpp:34
Evaluation cos(const Evaluation &value)
Definition: MathToolbox.hpp:383
unsigned invertQuadraticPolynomial(SolContainer &sol, Scalar a, Scalar b, Scalar c)
Invert a quadratic polynomial analytically.
Definition: PolynomialUtils.hpp:79
Evaluation acosh(const Evaluation &value)
Definition: MathToolbox.hpp:395
unsigned invertLinearPolynomial(SolContainer &sol, Scalar a, Scalar b)
Invert a linear polynomial analytically.
Definition: PolynomialUtils.hpp:51
Evaluation sqrt(const Evaluation &value)
Definition: MathToolbox.hpp:399
Evaluation asinh(const Evaluation &value)
Definition: MathToolbox.hpp:379
unsigned cubicRoots(SolContainer *sol, Scalar a, Scalar b, Scalar c, Scalar d)
Invert a cubic polynomial analytically.
Definition: PolynomialUtils.hpp:331
ReturnEval_< Evaluation1, Evaluation2 >::type atan2(const Evaluation1 &value1, const Evaluation2 &value2)
Definition: MathToolbox.hpp:363
auto scalarValue(const Evaluation &val) -> decltype(MathToolbox< Evaluation >::scalarValue(val))
Definition: MathToolbox.hpp:335
Evaluation acos(const Evaluation &value)
Definition: MathToolbox.hpp:387
Evaluation abs(const Evaluation &value)
Definition: MathToolbox.hpp:350
unsigned invertCubicPolynomial(SolContainer *sol, Scalar a, Scalar b, Scalar c, Scalar d)
Invert a cubic polynomial analytically.
Definition: PolynomialUtils.hpp:148
Evaluation cosh(const Evaluation &value)
Definition: MathToolbox.hpp:391
Evaluation sinh(const Evaluation &value)
Definition: MathToolbox.hpp:375
ReturnEval_< Evaluation1, Evaluation2 >::type pow(const Evaluation1 &base, const Evaluation2 &exp)
Definition: MathToolbox.hpp:416