37 const scalar
a = this->
a();
38 const scalar
b = this->
b();
39 const scalar
c = this->
c();
40 const scalar
d = this->
d();
44 <<
"Coefficients of the characteristic cubic polynomial:" <<
nl
60 const scalar
p = -(fma(-
a,
c, w) + fma(
b,
b/3.0, -w));
61 const scalar q =
b*
b*
b*scalar(2)/27 -
b*
c*
a/3 +
d*
a*
a;
62 const scalar numDiscr =
p*
p*
p/27 + q*q/4;
63 const scalar discr = (
mag(numDiscr) > VSMALL) ? numDiscr : 0;
66 const bool threeReal = discr < 0;
67 const bool oneRealTwoComplex = discr > 0;
69 (
mag(
p) >
sqrt(SMALL)) && !(threeReal || oneRealTwoComplex);
74 <<
"Numerical discriminant:" <<
tab << numDiscr <<
nl
75 <<
"Adjusted discriminant:" <<
tab << discr <<
nl
76 <<
"Number and types of the roots:" <<
nl
77 <<
"threeReal = " << threeReal <<
nl
78 <<
"oneRealTwoComplex = " << oneRealTwoComplex <<
nl
79 <<
"twoReal = " << twoReal <<
nl
80 <<
"oneReal = " << !(threeReal || oneRealTwoComplex || twoReal) <<
nl
84 static const scalar sqrt3 =
sqrt(3.0);
90 const scalar wCbRe = -q/2;
91 const scalar wCbIm =
sqrt(-discr);
92 const scalar wAbs =
cbrt(
hypot(wCbRe, wCbIm));
93 const scalar wArg =
atan2(wCbIm, wCbRe)/3;
94 const scalar wRe = wAbs*
cos(wArg);
95 const scalar wIm = wAbs*
sin(wArg);
99 x = -wRe -
mag(wIm)*sqrt3 -
b/3;
106 else if (oneRealTwoComplex)
108 const scalar wCb = -q/2 -
sign(q)*
sqrt(discr);
109 const scalar w =
cbrt(wCb);
110 const scalar t = w -
p/(3*w);
118 const scalar xRe = -t/2 -
b/3;
119 const scalar xIm = sqrt3/2*(w +
p/3/w);
155 <<
"#######" <<
endl;
Templated storage for the roots of polynomial equations, plus flags to indicate the nature of the roo...
void type(const direction i, const roots::type t)
Set the type of the i-th root.
Roots< 3 > roots() const
Return the roots of the cubic equation with no particular order.
Container to encapsulate various operations for linear equation of the forms with real coefficients:
Container to encapsulate various operations for quadratic equation of the forms with real coefficient...
dimensionedScalar sign(const dimensionedScalar &ds)
dimensionedScalar sin(const dimensionedScalar &ds)
messageStream Info
Information stream (stdout output on master, null elsewhere)
dimensionedScalar atan2(const dimensionedScalar &x, const dimensionedScalar &y)
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
dimensionedScalar hypot(const dimensionedScalar &x, const dimensionedScalar &y)
dimensionedScalar cbrt(const dimensionedScalar &ds)
dimensionedScalar cos(const dimensionedScalar &ds)
constexpr char nl
The newline '\n' character (0x0a)
constexpr char tab
The tab '\t' character(0x09)