Commit 5e7a7661 authored by Will Bainbridge's avatar Will Bainbridge Committed by Andrew Heather
Browse files

BUG: cubicEqn, quadraticEqn: Correction to repeated roots

Also extended the cubic equation test routine and modified the error
methods so that they more accurately generate the round of error of
evaluation.

This resolves bug report https://bugs.openfoam.org/view.php?id=3015
parent d99707fa
......@@ -51,7 +51,7 @@ void test(const Type& polynomialEqn, const scalar tol)
case roots::real:
v[i] = polynomialEqn.value(r[i]);
e[i] = polynomialEqn.error(r[i]);
ok = ok && mag(v[i]) < tol*mag(e[i]);
ok = ok && mag(v[i]) <= tol*mag(e[i]);
break;
case roots::posInf:
v[i] = + inf;
......@@ -79,8 +79,10 @@ void test(const Type& polynomialEqn, const scalar tol)
int main()
{
const int nTests = 1000000;
for (int t = 0; t < nTests; ++ t)
const scalar tol = 5;
const label nTests = 1000000;
for (label t = 0; t < nTests; ++ t)
{
test
(
......@@ -91,11 +93,26 @@ int main()
randomScalar(1e-50, 1e+50),
randomScalar(1e-50, 1e+50)
),
100
tol
);
}
Info << nTests << " random cubics tested" << endl;
Info << nTests << " cubics tested" << endl;
const label coeffMin = -9, coeffMax = 10, nCoeff = coeffMax - coeffMin;
for (label a = coeffMin; a < coeffMax; ++ a)
{
for (label b = coeffMin; b < coeffMax; ++ b)
{
for (label c = coeffMin; c < coeffMax; ++ c)
{
for (label d = coeffMin; d < coeffMax; ++ d)
{
test(cubicEqn(a, b, c, d), tol);
}
}
}
}
Info<< nCoeff*nCoeff*nCoeff*nCoeff << " integer cubics tested" << endl;
return 0;
}
......@@ -81,7 +81,7 @@ Foam::Roots<3> Foam::cubicEqn::roots() const
// This is assumed not to over- or under-flow. If it does, all bets are off.
const scalar p = c*a - b*b/3;
const scalar q = b*b*b*(2.0/27.0) - b*c*a/3 + d*a*a;
const scalar q = b*b*b*scalar(2)/27 - b*c*a/3 + d*a*a;
const scalar disc = p*p*p/27 + q*q/4;
// How many roots of what types are available?
......@@ -96,7 +96,7 @@ Foam::Roots<3> Foam::cubicEqn::roots() const
if (oneReal)
{
const Roots<1> r = linearEqn(- a, b/3).roots();
const Roots<1> r = linearEqn(a, b/3).roots();
return Roots<3>(r.type(0), r[0]);
}
else if (twoReal)
......
......@@ -114,7 +114,10 @@ inline Foam::scalar Foam::cubicEqn::derivative(const scalar x) const
inline Foam::scalar Foam::cubicEqn::error(const scalar x) const
{
return mag(SMALL*x*(x*(x*3*a() + 2*b()) + c()));
return
SMALL*magSqr(x)*(mag(x*a()) + mag(b()))
+ SMALL*mag(x)*(mag(x*(x*a() + b())) + mag(c()))
+ SMALL*(mag(x*(x*(x*a() + b()) + c())) + mag(d()));
}
......
......@@ -82,7 +82,7 @@ inline Foam::scalar Foam::linearEqn::derivative(const scalar x) const
inline Foam::scalar Foam::linearEqn::error(const scalar x) const
{
return mag(SMALL*x*a());
return SMALL*(mag(x*a()) + mag(b()));
}
......
......@@ -71,7 +71,7 @@ Foam::Roots<2> Foam::quadraticEqn::roots() const
if (oneReal)
{
const Roots<1> r = linearEqn(- a, b/2).roots();
const Roots<1> r = linearEqn(a, b/2).roots();
return Roots<2>(r, r);
}
else if (twoReal)
......
......@@ -100,7 +100,9 @@ inline Foam::scalar Foam::quadraticEqn::derivative(const scalar x) const
inline Foam::scalar Foam::quadraticEqn::error(const scalar x) const
{
return mag(SMALL*x*(x*2*a() + b()));
return
SMALL*mag(x)*(mag(x*a()) + mag(b()))
+ SMALL*(mag(x*(x*a() + b())) + mag(c()));
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment