diff --git a/applications/test/cubicEqn/Test-cubicEqn.C b/applications/test/cubicEqn/Test-cubicEqn.C
deleted file mode 100644
index aa0579a41940d23aac180a523d7b871a73f0d84f..0000000000000000000000000000000000000000
--- a/applications/test/cubicEqn/Test-cubicEqn.C
+++ /dev/null
@@ -1,118 +0,0 @@
-#include <ctime>
-#include <random>
-
-#include "cubicEqn.H"
-#include "IOstreams.H"
-#include "stringList.H"
-
-using namespace Foam;
-
-scalar randomScalar(const scalar min, const scalar max)
-{
- static_assert
- (
- sizeof(long) == sizeof(scalar),
- "Scalar and long are not the same size"
- );
- static std::default_random_engine generator(std::time(0));
- static std::uniform_int_distribution<long>
- distribution
- (
- std::numeric_limits<long>::min(),
- std::numeric_limits<long>::max()
- );
- scalar x;
- do
- {
- long i = distribution(generator);
- x = reinterpret_cast<scalar&>(i);
- }
- while (min > mag(x) || mag(x) > max || !std::isfinite(x));
- return x;
-};
-
-template <class Type>
-void test(const Type& polynomialEqn, const scalar tol)
-{
- Roots<Type::nComponents - 1> r = polynomialEqn.roots();
-
- const scalar nan = std::numeric_limits<scalar>::quiet_NaN();
- const scalar inf = std::numeric_limits<scalar>::infinity();
-
- FixedList<label, Type::nComponents - 1> t;
- FixedList<scalar, Type::nComponents - 1> v(nan);
- FixedList<scalar, Type::nComponents - 1> e(nan);
- bool ok = true;
- forAll(r, i)
- {
- t[i] = r.type(i);
- switch (t[i])
- {
- case roots::real:
- v[i] = polynomialEqn.value(r[i]);
- e[i] = polynomialEqn.error(r[i]);
- ok = ok && mag(v[i]) <= tol*mag(e[i]);
- break;
- case roots::posInf:
- v[i] = + inf;
- e[i] = nan;
- break;
- case roots::negInf:
- v[i] = - inf;
- e[i] = nan;
- break;
- default:
- v[i] = e[i] = nan;
- break;
- }
- }
-
- if (!ok)
- {
- Info<< "Coeffs: " << polynomialEqn << endl
- << " Types: " << t << endl
- << " Roots: " << r << endl
- << "Values: " << v << endl
- << "Errors: " << e << endl << endl;
- }
-}
-
-int main()
-{
- const scalar tol = 5;
-
- const label nTests = 1000000;
- for (label t = 0; t < nTests; ++ t)
- {
- test
- (
- cubicEqn
- (
- randomScalar(1e-50, 1e+50),
- randomScalar(1e-50, 1e+50),
- randomScalar(1e-50, 1e+50),
- randomScalar(1e-50, 1e+50)
- ),
- tol
- );
- }
- Info << nTests << " random 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;
-}
diff --git a/applications/test/cubicEqn/Make/files b/applications/test/polynomialEqns/cubicEqn/Make/files
similarity index 100%
rename from applications/test/cubicEqn/Make/files
rename to applications/test/polynomialEqns/cubicEqn/Make/files
diff --git a/applications/test/cubicEqn/Make/options b/applications/test/polynomialEqns/cubicEqn/Make/options
similarity index 100%
rename from applications/test/cubicEqn/Make/options
rename to applications/test/polynomialEqns/cubicEqn/Make/options
diff --git a/applications/test/polynomialEqns/cubicEqn/Test-cubicEqn.C b/applications/test/polynomialEqns/cubicEqn/Test-cubicEqn.C
new file mode 100644
index 0000000000000000000000000000000000000000..df510df1959c9fd5eae9eb759297e0d1a2db82b7
--- /dev/null
+++ b/applications/test/polynomialEqns/cubicEqn/Test-cubicEqn.C
@@ -0,0 +1,339 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2018 OpenFOAM Foundation
+ Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
+
+Application
+ Test-cubicEqn
+
+Description
+ Tests for \c cubicEqn constructors, and member functions.
+
+\*---------------------------------------------------------------------------*/
+
+#include <ctime>
+#include <random>
+
+#include "cubicEqn.H"
+#include "vector.H"
+#include "scalarList.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Total number of unit tests
+unsigned nTest_ = 0;
+
+
+// Total number of failed unit tests
+unsigned nFail_ = 0;
+
+
+// Compare two floating point types, and print output.
+// Do ++nFail_ if values of two objects are not equal within a given tolerance.
+// The function is converted from PEP-485.
+void cmp
+(
+ const word& msg,
+ const scalar& x,
+ const scalar& y,
+ const scalar relTol = 1e-8, //<! are values the same within 8 decimals
+ const scalar absTol = 0 //<! useful for cmps near zero
+)
+{
+ Info<< msg << x << endl;
+
+ unsigned nFail = 0;
+
+ if (max(absTol, relTol*max(mag(x), mag(y))) < mag(x - y))
+ {
+ ++nFail;
+ }
+
+ if (nFail)
+ {
+ Info<< nl
+ << " #### Fail in " << nFail << " comps ####" << nl << endl;
+ ++nFail_;
+ }
+ ++nTest_;
+}
+
+
+// Create each constructor of cubicEqn, and print output
+void test_constructors()
+{
+ #if 0
+ {
+ Info<< "# Construct null:" << nl;
+ const cubicEqn T();
+ Info<< T << endl;
+ }
+ #endif
+
+ {
+ Info<< "# Construct initialized to zero:" << nl;
+ const cubicEqn T(Zero);
+ Info<< T << endl;
+ }
+ {
+ Info<< "# Construct from components:" << nl;
+ const cubicEqn T(1, 2, 3, 4);
+ Info<< T << endl;
+ }
+}
+
+
+// Execute each member function of cubicEqn, and print output
+void test_member_funcs()
+{
+ const cubicEqn T(1, -6, -2, 12);
+ const scalar x = 9.7;
+
+ Info<< "# Operands: " << nl
+ << " cubicEqn = " << T << nl
+ << " x = " << x << endl;
+
+
+ {
+ Info<< "# Access:" << nl;
+ cmp(" a = ", 1, T.a());
+ cmp(" b = ", -6, T.b());
+ cmp(" c = ", -2, T.c());
+ cmp(" d = ", 12, T.d());
+ }
+ {
+ Info<< "# Evaluate:" << nl;
+ cmp(" T.value(x) = ", T.value(x), 340.73299999999983);
+ cmp(" T.derivative(x) = ", T.derivative(x), 163.87);
+ cmp
+ (
+ " T.error(x) = ",
+ T.error(x),
+ 2.1854789999999994e-13,
+ 1e-8,
+ 1e-8
+ );
+
+ const vector rts(6, 1.41421356, -1.41421356);
+ forAll(rts, i)
+ {
+ cmp(" T.roots() = ", T.roots()[i], rts[i]);
+ }
+ }
+
+ {
+ Info<< "# Special cases for the root evaluation: " << nl;
+
+ {
+ Info<< " Three distinct real roots:" << nl;
+ const cubicEqn Tc(1.0, -4.0, 1.0, 6.0);
+ const vector rts(3, 2, -1);
+ forAll(rts, i)
+ {
+ cmp(" root = ", Tc.roots()[i], rts[i]);
+ }
+ }
+ {
+ Info<< " One real root and one complex conjugate-pair root:" << nl;
+ const cubicEqn Tc(1.0, 2.0, 3.0, 4.0);
+ const vector rts(-1.65062919, -0.1746854, 1.54686889);
+ forAll(rts, i)
+ {
+ cmp(" root = ", Tc.roots()[i], rts[i], 1e-8, 1e-8);
+ }
+ }
+ {
+ Info<< " Two identical real roots and a distinct real root:" << nl;
+ const cubicEqn Tc(1.0, -5.0, 8.0, -4.0);
+ const vector rts(2, 2, 1);
+ forAll(rts, i)
+ {
+ cmp(" root = ", Tc.roots()[i], rts[i]);
+ }
+ }
+ {
+ Info<< " Three identical real roots:" << nl;
+ const cubicEqn Tc(1.0, -6.0, 12.0, -8.0);
+ const vector rts(2, 2, 2);
+ forAll(rts, i)
+ {
+ cmp(" root = ", Tc.roots()[i], rts[i]);
+ }
+ }
+ {
+ Info<< " Arbitrary roots with non-unity leading term:" << nl;
+ const cubicEqn Tc(-8.59, 6.77, -0.2, 8.0);
+ const vector rts(1.31167947, -0.26177687, 0.80093099);
+ forAll(rts, i)
+ {
+ cmp(" root = ", Tc.roots()[i], rts[i]);
+ }
+ }
+ }
+}
+
+
+scalar randomScalar(const scalar min, const scalar max)
+{
+ static_assert
+ (
+ sizeof(long) == sizeof(scalar),
+ "Scalar and long are not the same size"
+ );
+ static std::default_random_engine generator(std::time(0));
+ static std::uniform_int_distribution<long>
+ distribution
+ (
+ std::numeric_limits<long>::min(),
+ std::numeric_limits<long>::max()
+ );
+
+ scalar x = 0;
+
+ do
+ {
+ x = scalar(distribution(generator));
+ }
+ while (min > mag(x) || mag(x) > max || !std::isfinite(x));
+
+ return x;
+};
+
+
+template <class Type>
+void test(const Type& polynomialEqn, const scalar tol)
+{
+ Roots<Type::nComponents - 1> r = polynomialEqn.roots();
+
+ const scalar nan = std::numeric_limits<scalar>::quiet_NaN();
+ const scalar inf = std::numeric_limits<scalar>::infinity();
+
+ FixedList<label, Type::nComponents - 1> t;
+ FixedList<scalar, Type::nComponents - 1> v(nan);
+ FixedList<scalar, Type::nComponents - 1> e(nan);
+ bool ok = true;
+ forAll(r, i)
+ {
+ t[i] = r.type(i);
+ switch (t[i])
+ {
+ case roots::real:
+ v[i] = polynomialEqn.value(r[i]);
+ e[i] = polynomialEqn.error(r[i]);
+ ok = ok && mag(v[i]) <= tol*mag(e[i]);
+ break;
+ case roots::posInf:
+ v[i] = + inf;
+ e[i] = nan;
+ break;
+ case roots::negInf:
+ v[i] = - inf;
+ e[i] = nan;
+ break;
+ default:
+ v[i] = e[i] = nan;
+ break;
+ }
+ }
+
+ if (!ok)
+ {
+ Info<< "Coeffs: " << polynomialEqn << endl
+ << " Types: " << t << endl
+ << " Roots: " << r << endl
+ << "Values: " << v << endl
+ << "Errors: " << e << endl << endl;
+ ++nFail_;
+ }
+}
+
+
+// * * * * * * * * * * * * * * * Main Program * * * * * * * * * * * * * * * //
+
+int main()
+{
+ Info<< nl << " ## Test constructors: ##" << nl;
+ test_constructors();
+
+ Info<< nl << " ## Test member functions: ##" << nl;
+ test_member_funcs();
+
+
+ // Pre-v2006 tests
+ const scalar tol = 5;
+
+ const label nTests = 1000000;
+ for (label t = 0; t < nTests; ++ t)
+ {
+ test
+ (
+ cubicEqn
+ (
+ randomScalar(1e-50, 1e+50),
+ randomScalar(1e-50, 1e+50),
+ randomScalar(1e-50, 1e+50),
+ randomScalar(1e-50, 1e+50)
+ ),
+ tol
+ );
+ ++nTest_;
+ }
+ Info << nTests << " scalar cubic equations were 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);
+ ++nTest_;
+ }
+ }
+ }
+ }
+ Info<< nCoeff*nCoeff*nCoeff*nCoeff
+ << " label cubic equations were tested." << endl;
+
+
+ if (nFail_)
+ {
+ Info<< nl << " #### "
+ << "Failed in " << nFail_ << " tests "
+ << "out of total " << nTest_ << " tests "
+ << "####\n" << endl;
+ return 1;
+ }
+
+ Info<< nl << " #### Passed all " << nTest_ <<" tests ####\n" << endl;
+ return 0;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/test/polynomialEqns/linearEqn/Make/files b/applications/test/polynomialEqns/linearEqn/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..d897b48a2e87aed3750ffe90be5fb0320395c631
--- /dev/null
+++ b/applications/test/polynomialEqns/linearEqn/Make/files
@@ -0,0 +1,3 @@
+Test-linearEqn.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-linearEqn
diff --git a/applications/test/polynomialEqns/linearEqn/Make/options b/applications/test/polynomialEqns/linearEqn/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..4c3dd783cb4170feefb3f5385510a83257b43b18
--- /dev/null
+++ b/applications/test/polynomialEqns/linearEqn/Make/options
@@ -0,0 +1,3 @@
+EXE_INC =
+
+EXE_LIBS =
diff --git a/applications/test/polynomialEqns/linearEqn/Test-linearEqn.C b/applications/test/polynomialEqns/linearEqn/Test-linearEqn.C
new file mode 100644
index 0000000000000000000000000000000000000000..683a1b161d9c4fee33ed203634c6eb572d23cc97
--- /dev/null
+++ b/applications/test/polynomialEqns/linearEqn/Test-linearEqn.C
@@ -0,0 +1,270 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
+
+Application
+ Test-linearEqn
+
+Description
+ Tests for \c linearEqn constructors, and member functions.
+
+\*---------------------------------------------------------------------------*/
+
+#include <ctime>
+#include <random>
+
+#include "linearEqn.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Total number of unit tests
+unsigned nTest_ = 0;
+
+
+// Total number of failed unit tests
+unsigned nFail_ = 0;
+
+
+// Compare two floating point types, and print output.
+// Do ++nFail_ if values of two objects are not equal within a given tolerance.
+// The function is converted from PEP-485.
+void cmp
+(
+ const word& msg,
+ const scalar& x,
+ const scalar& y,
+ const scalar relTol = 1e-8, //<! are values the same within 8 decimals
+ const scalar absTol = 0 //<! useful for cmps near zero
+)
+{
+ Info<< msg << x << endl;
+
+ unsigned nFail = 0;
+
+ if (max(absTol, relTol*max(mag(x), mag(y))) < mag(x - y))
+ {
+ ++nFail;
+ }
+
+ if (nFail)
+ {
+ Info<< nl
+ << " #### Fail in " << nFail << " comps ####" << nl << endl;
+ ++nFail_;
+ }
+ ++nTest_;
+}
+
+
+// Create each constructor of linearEqn, and print output
+void test_constructors()
+{
+ #if 0
+ {
+ Info<< "# Construct null:" << nl;
+ const linearEqn T();
+ Info<< T << endl;
+ }
+ #endif
+
+ {
+ Info<< "# Construct initialized to zero:" << nl;
+ const linearEqn T(Zero);
+ Info<< T << endl;
+ }
+ {
+ Info<< "# Construct from components:" << nl;
+ const linearEqn T(1, 2);
+ Info<< T << endl;
+ }
+}
+
+
+// Execute each member function of linearEqn, and print output
+void test_member_funcs()
+{
+ const linearEqn T(1, -6);
+ const scalar x = 9.7;
+
+ Info<< "# Operands: " << nl
+ << " linearEqn = " << T << nl
+ << " x = " << x << endl;
+
+
+ {
+ Info<< "# Access:" << nl;
+ cmp(" a = ", 1, T.a());
+ cmp(" b = ", -6, T.b());
+ }
+ {
+ Info<< "# Evaluate:" << nl;
+ cmp(" T.value(x) = ", T.value(x), 3.6999999999999993);
+ cmp(" T.derivative(x) = ", T.derivative(x), 1);
+ cmp
+ (
+ " T.error(x) = ",
+ T.error(x),
+ 1.5699999999999999e-15,
+ 1e-8,
+ 1e-8
+ );
+ cmp(" T.roots() = ", T.roots()[0], 6);
+ }
+}
+
+
+scalar randomScalar(const scalar min, const scalar max)
+{
+ static_assert
+ (
+ sizeof(long) == sizeof(scalar),
+ "Scalar and long are not the same size"
+ );
+ static std::default_random_engine generator(std::time(0));
+ static std::uniform_int_distribution<long>
+ distribution
+ (
+ std::numeric_limits<long>::min(),
+ std::numeric_limits<long>::max()
+ );
+
+ scalar x = 0;
+
+ do
+ {
+ x = scalar(distribution(generator));
+ }
+ while (min > mag(x) || mag(x) > max || !std::isfinite(x));
+
+ return x;
+};
+
+
+template <class Type>
+void test(const Type& polynomialEqn, const scalar tol)
+{
+ Roots<Type::nComponents - 1> r = polynomialEqn.roots();
+
+ const scalar nan = std::numeric_limits<scalar>::quiet_NaN();
+ const scalar inf = std::numeric_limits<scalar>::infinity();
+
+ FixedList<label, Type::nComponents - 1> t;
+ FixedList<scalar, Type::nComponents - 1> v(nan);
+ FixedList<scalar, Type::nComponents - 1> e(nan);
+ bool ok = true;
+ forAll(r, i)
+ {
+ t[i] = r.type(i);
+ switch (t[i])
+ {
+ case roots::real:
+ v[i] = polynomialEqn.value(r[i]);
+ e[i] = polynomialEqn.error(r[i]);
+ ok = ok && mag(v[i]) <= tol*mag(e[i]);
+ break;
+ case roots::posInf:
+ v[i] = + inf;
+ e[i] = nan;
+ break;
+ case roots::negInf:
+ v[i] = - inf;
+ e[i] = nan;
+ break;
+ default:
+ v[i] = e[i] = nan;
+ break;
+ }
+ }
+
+ if (!ok)
+ {
+ Info<< "Coeffs: " << polynomialEqn << endl
+ << " Types: " << t << endl
+ << " Roots: " << r << endl
+ << "Values: " << v << endl
+ << "Errors: " << e << endl << endl;
+ ++nFail_;
+ }
+}
+
+
+// * * * * * * * * * * * * * * * Main Program * * * * * * * * * * * * * * * //
+
+int main()
+{
+ Info<< nl << " ## Test constructors: ##" << nl;
+ test_constructors();
+
+ Info<< nl << " ## Test member functions: ##" << nl;
+ test_member_funcs();
+
+
+ // Sanity checks
+ const scalar tol = 5;
+
+ const label nTests = 1000000;
+ for (label t = 0; t < nTests; ++t)
+ {
+ test
+ (
+ linearEqn
+ (
+ randomScalar(1e-50, 1e+50),
+ randomScalar(1e-50, 1e+50)
+ ),
+ tol
+ );
+ ++nTest_;
+ }
+ Info << nTests << " scalar linear equations were 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)
+ {
+ test(linearEqn(a, b), tol);
+ ++nTest_;
+ }
+ }
+ Info<< nCoeff*nCoeff << " label linear equations were tested." << endl;
+
+
+ if (nFail_)
+ {
+ Info<< nl << " #### "
+ << "Failed in " << nFail_ << " tests "
+ << "out of total " << nTest_ << " tests "
+ << "####\n" << endl;
+ return 1;
+ }
+
+ Info<< nl << " #### Passed all " << nTest_ <<" tests ####\n" << endl;
+ return 0;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/test/polynomialEqns/quadraticEqn/Make/files b/applications/test/polynomialEqns/quadraticEqn/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..997b2f63cb518a617da86f4a8b3f320ff4da615d
--- /dev/null
+++ b/applications/test/polynomialEqns/quadraticEqn/Make/files
@@ -0,0 +1,3 @@
+Test-quadraticEqn.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-quadraticEqn
diff --git a/applications/test/polynomialEqns/quadraticEqn/Make/options b/applications/test/polynomialEqns/quadraticEqn/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..4c3dd783cb4170feefb3f5385510a83257b43b18
--- /dev/null
+++ b/applications/test/polynomialEqns/quadraticEqn/Make/options
@@ -0,0 +1,3 @@
+EXE_INC =
+
+EXE_LIBS =
diff --git a/applications/test/polynomialEqns/quadraticEqn/Test-quadraticEqn.C b/applications/test/polynomialEqns/quadraticEqn/Test-quadraticEqn.C
new file mode 100644
index 0000000000000000000000000000000000000000..61e28e394ce49b0f472b07a3ec51f895d7d06134
--- /dev/null
+++ b/applications/test/polynomialEqns/quadraticEqn/Test-quadraticEqn.C
@@ -0,0 +1,321 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
+
+Application
+ Test-quadraticEqn
+
+Description
+ Tests for \c quadraticEqn constructors, and member functions.
+
+\*---------------------------------------------------------------------------*/
+
+#include <ctime>
+#include <random>
+
+#include "quadraticEqn.H"
+#include "vector2D.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Total number of unit tests
+unsigned nTest_ = 0;
+
+
+// Total number of failed unit tests
+unsigned nFail_ = 0;
+
+
+// Compare two floating point types, and print output.
+// Do ++nFail_ if values of two objects are not equal within a given tolerance.
+// The function is converted from PEP-485.
+void cmp
+(
+ const word& msg,
+ const scalar& x,
+ const scalar& y,
+ const scalar relTol = 1e-8, //<! are values the same within 8 decimals
+ const scalar absTol = 0 //<! useful for cmps near zero
+)
+{
+ Info<< msg << x << endl;
+
+ unsigned nFail = 0;
+
+ if (max(absTol, relTol*max(mag(x), mag(y))) < mag(x - y))
+ {
+ ++nFail;
+ }
+
+ if (nFail)
+ {
+ Info<< nl
+ << " #### Fail in " << nFail << " comps ####" << nl << endl;
+ ++nFail_;
+ }
+ ++nTest_;
+}
+
+
+// Create each constructor of quadraticEqn, and print output
+void test_constructors()
+{
+ #if 0
+ {
+ Info<< "# Construct null:" << nl;
+ const quadraticEqn T();
+ Info<< T << endl;
+ }
+ #endif
+
+ {
+ Info<< "# Construct initialized to zero:" << nl;
+ const quadraticEqn T(Zero);
+ Info<< T << endl;
+ }
+ {
+ Info<< "# Construct from components:" << nl;
+ const quadraticEqn T(1, 2, 3);
+ Info<< T << endl;
+ }
+}
+
+
+// Execute each member function of quadraticEqn, and print output
+void test_member_funcs()
+{
+ const quadraticEqn T(1, -6, -2);
+ const scalar x = 9.7;
+
+ Info<< "# Operands: " << nl
+ << " quadraticEqn = " << T << nl
+ << " x = " << x << endl;
+
+ {
+ Info<< "# Access:" << nl;
+ cmp(" a = ", 1, T.a());
+ cmp(" b = ", -6, T.b());
+ cmp(" c = ", -2, T.c());
+ }
+ {
+ Info<< "# Evaluate:" << nl;
+ cmp(" T.value(x) = ", T.value(x), 33.88999999999999);
+ cmp(" T.derivative(x) = ", T.derivative(x), 13.399999999999999);
+ cmp
+ (
+ " T.error(x) = ",
+ T.error(x),
+ 1.9017999999999998e-14,
+ 1e-8,
+ 1e-8
+ );
+
+ const vector2D rts(6.31662479, -0.31662479);
+ forAll(rts, i)
+ {
+ cmp(" T.roots() = ", T.roots()[i], rts[i]);
+ }
+ }
+ {
+ Info<< "# Special cases for the root evaluation: " << nl;
+
+ {
+ Info<< " Two distinct real roots:" << nl;
+ const quadraticEqn Tc(1.0, -11.0, 24.0);
+ const vector2D rts(8, 3);
+ forAll(rts, i)
+ {
+ cmp(" root = ", Tc.roots()[i], rts[i]);
+ }
+ }
+ {
+ Info<< " Two identical real roots:" << nl;
+ const quadraticEqn Tc(1.0, -8.0, 16.0);
+ const vector2D rts(4, 4);
+ forAll(rts, i)
+ {
+ cmp(" root = ", Tc.roots()[i], rts[i]);
+ }
+ }
+ {
+ Info<< " One complex conjugate-pair root:" << nl;
+ const quadraticEqn Tc(2.0, -2.0, 1.0);
+ const vector2D rts(0.5, -0.5);
+ forAll(rts, i)
+ {
+ cmp(" root = ", Tc.roots()[i], rts[i], 1e-8, 1e-8);
+ }
+ }
+ {
+ Info<< " Cancellation problem:" << nl;
+ const quadraticEqn Tc(1.0, 68.5, 0.1);
+ const vector2D rts(-6.84985401e+01, -1.45988513e-03);
+ forAll(rts, i)
+ {
+ cmp(" root = ", Tc.roots()[i], rts[i], 1e-8, 1e-8);
+ }
+ }
+ }
+}
+
+
+scalar randomScalar(const scalar min, const scalar max)
+{
+ static_assert
+ (
+ sizeof(long) == sizeof(scalar),
+ "Scalar and long are not the same size"
+ );
+ static std::default_random_engine generator(std::time(0));
+ static std::uniform_int_distribution<long>
+ distribution
+ (
+ std::numeric_limits<long>::min(),
+ std::numeric_limits<long>::max()
+ );
+
+ scalar x = 0;
+
+ do
+ {
+ x = scalar(distribution(generator));
+ }
+ while (min > mag(x) || mag(x) > max || !std::isfinite(x));
+
+ return x;
+};
+
+
+template <class Type>
+void test(const Type& polynomialEqn, const scalar tol)
+{
+ Roots<Type::nComponents - 1> r = polynomialEqn.roots();
+
+ const scalar nan = std::numeric_limits<scalar>::quiet_NaN();
+ const scalar inf = std::numeric_limits<scalar>::infinity();
+
+ FixedList<label, Type::nComponents - 1> t;
+ FixedList<scalar, Type::nComponents - 1> v(nan);
+ FixedList<scalar, Type::nComponents - 1> e(nan);
+ bool ok = true;
+ forAll(r, i)
+ {
+ t[i] = r.type(i);
+ switch (t[i])
+ {
+ case roots::real:
+ v[i] = polynomialEqn.value(r[i]);
+ e[i] = polynomialEqn.error(r[i]);
+ ok = ok && mag(v[i]) <= tol*mag(e[i]);
+ break;
+ case roots::posInf:
+ v[i] = + inf;
+ e[i] = nan;
+ break;
+ case roots::negInf:
+ v[i] = - inf;
+ e[i] = nan;
+ break;
+ default:
+ v[i] = e[i] = nan;
+ break;
+ }
+ }
+
+ if (!ok)
+ {
+ Info<< "Coeffs: " << polynomialEqn << endl
+ << " Types: " << t << endl
+ << " Roots: " << r << endl
+ << "Values: " << v << endl
+ << "Errors: " << e << endl << endl;
+ ++nFail_;
+ }
+}
+
+
+// * * * * * * * * * * * * * * * Main Program * * * * * * * * * * * * * * * //
+
+int main()
+{
+ Info<< nl << " ## Test constructors: ##" << nl;
+ test_constructors();
+
+ Info<< nl << " ## Test member functions: ##" << nl;
+ test_member_funcs();
+
+
+ // Sanity checks
+ const scalar tol = 5;
+
+ const label nTests = 1000000;
+ for (label t = 0; t < nTests; ++t)
+ {
+ test
+ (
+ quadraticEqn
+ (
+ randomScalar(1e-50, 1e+50),
+ randomScalar(1e-50, 1e+50),
+ randomScalar(1e-50, 1e+50)
+ ),
+ tol
+ );
+ ++nTest_;
+ }
+ Info << nTests << " scalar quadratic equations were 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)
+ {
+ test(quadraticEqn(a, b, c), tol);
+ ++nTest_;
+ }
+ }
+ }
+ Info<< nCoeff*nCoeff*nCoeff
+ << " label quadratic equations were tested." << endl;
+
+
+ if (nFail_)
+ {
+ Info<< nl << " #### "
+ << "Failed in " << nFail_ << " tests "
+ << "out of total " << nTest_ << " tests "
+ << "####\n" << endl;
+ return 1;
+ }
+
+ Info<< nl << " #### Passed all " << nTest_ <<" tests ####\n" << endl;
+ return 0;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.C b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.C
index 35ad8934417d6bb94640853faa772fcc13c82dcc..fc69445e9708b152c56ff12bb530436375008004 100644
--- a/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.C
+++ b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.C
@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2017 OpenFOAM Foundation
+ Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -33,132 +34,104 @@ License
Foam::Roots<3> Foam::cubicEqn::roots() const
{
- /*
-
- This function solves a cubic equation of the following form:
-
- a*x^3 + b*x^2 + c*x + d = 0
- x^3 + B*x^2 + C*x + D = 0
-
- The following two substitutions are used:
-
- x = t - B/3
- t = w - P/3/w
-
- This reduces the problem to a quadratic in w^3.
-
- w^6 + Q*w^3 - P = 0
-
- Where Q and P are given in the code below.
-
- The properties of the cubic can be related to the properties of this
- quadratic in w^3. If it has a repeated root a zero, the cubic has a tripl
- root. If it has a repeated root not at zero, the cubic has two real roots,
- one repeated and one not. If it has two complex roots, the cubic has three
- real roots. If it has two real roots, then the cubic has one real root and
- two complex roots.
-
- This is solved for the most numerically accurate value of w^3. See the
- quadratic function for details on how to pick a value. This single value of
- w^3 can yield up to three cube roots for w, which relate to the three
- solutions for x.
-
- Only a single root, or pair of conjugate roots, is directly evaluated; the
- one, or ones with the lowest relative numerical error. Root identities are
- then used to recover the remaining roots, possibly utilising a quadratic
- and/or linear solution. This seems to be a good way of maintaining the
- accuracy of roots at very different magnitudes.
-
- */
-
const scalar a = this->a();
const scalar b = this->b();
const scalar c = this->c();
const scalar d = this->d();
- if (a == 0)
+ // Check the leading term in the cubic eqn exists
+ if (mag(a) < VSMALL)
{
return Roots<3>(quadraticEqn(b, c, d).roots(), roots::nan, 0);
}
- // This is assumed not to over- or under-flow. If it does, all bets are off.
- const scalar p = c*a - b*b/3;
+ // (JLM:p. 2246) [p = a*c - b*b/3]
+ const scalar w = a*c;
+ const scalar p = -(fma(-a, c, w) + fma(b, b/3.0, -w));
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;
+ const scalar numDiscr = p*p*p/27 + q*q/4;
+ const scalar discr = (mag(numDiscr) > SMALL) ? numDiscr : 0;
- // How many roots of what types are available?
- const bool oneReal = disc == 0 && p == 0;
- const bool twoReal = disc == 0 && p != 0;
- const bool threeReal = disc < 0;
- //const bool oneRealTwoComplex = disc > 0;
+ // Determine the number and types of the roots
+ const bool threeReal = discr < 0;
+ const bool oneRealTwoComplex = discr > 0;
+ const bool twoReal = //p != 0; && discr == 0;
+ (mag(p) > VSMALL) && !(threeReal || oneRealTwoComplex);
+ // const bool oneReal = p == 0 && discr == 0;
static const scalar sqrt3 = sqrt(3.0);
- scalar x;
+ scalar x = 0;
- if (oneReal)
- {
- const Roots<1> r = linearEqn(a, b/3).roots();
- return Roots<3>(r.type(0), r[0]);
- }
- else if (twoReal)
- {
- if (q*b > 0)
- {
- x = - 2*cbrt(q/2) - b/3;
- }
- else
- {
- x = cbrt(q/2) - b/3;
- const Roots<1> r = linearEqn(- a, x).roots();
- return Roots<3>(Roots<2>(r, r), linearEqn(x*x, a*d).roots());
- }
- }
- else if (threeReal)
+ if (threeReal)
{
- const scalar wCbRe = - q/2, wCbIm = sqrt(- disc);
+ const scalar wCbRe = -q/2;
+ const scalar wCbIm = sqrt(-discr);
const scalar wAbs = cbrt(hypot(wCbRe, wCbIm));
const scalar wArg = atan2(wCbIm, wCbRe)/3;
- const scalar wRe = wAbs*cos(wArg), wIm = wAbs*sin(wArg);
+ const scalar wRe = wAbs*cos(wArg);
+ const scalar wIm = wAbs*sin(wArg);
+
if (b > 0)
{
- x = - wRe - mag(wIm)*sqrt3 - b/3;
+ x = -wRe - mag(wIm)*sqrt3 - b/3;
}
else
{
x = 2*wRe - b/3;
}
}
- else // if (oneRealTwoComplex)
+ else if (oneRealTwoComplex)
{
- const scalar wCb = - q/2 - sign(q)*sqrt(disc);
+ const scalar wCb = -q/2 - sign(q)*sqrt(discr);
const scalar w = cbrt(wCb);
const scalar t = w - p/(3*w);
+
if (p + t*b < 0)
{
x = t - b/3;
}
else
{
- const scalar xRe = - t/2 - b/3, xIm = sqrt3/2*(w + p/3/w);
- x = - a*a*d/(xRe*xRe + xIm*xIm);
+ const scalar xRe = -t/2 - b/3;
+ const scalar xIm = sqrt3/2*(w + p/3/w);
- // This form of solving for the remaining roots seems more stable
- // for this case. This has been determined by trial and error.
return
Roots<3>
(
- linearEqn(- a, x).roots(),
- quadraticEqn(a*x, x*x + b*x, - a*d).roots()
+ Roots<1>(roots::real, -a*d/(xRe*xRe + xIm*xIm)),
+ Roots<2>
+ (
+ Roots<1>(roots::complex, xRe),
+ Roots<1>(roots::complex, xIm)
+ )
);
}
}
+ else if (twoReal)
+ {
+ if (q*b > 0)
+ {
+ x = -2*cbrt(q/2) - b/3;
+ }
+ else
+ {
+ x = cbrt(q/2) - b/3;
+ const Roots<1> r(linearEqn(-a, x).roots());
+ return Roots<3>(Roots<2>(r, r), linearEqn(x*x, a*d).roots());
+ }
+ }
+ else // (oneReal)
+ {
+ const Roots<1> r(linearEqn(a, b/3).roots());
+ return Roots<3>(r.type(0), r[0]);
+ }
return
Roots<3>
(
- linearEqn(- a, x).roots(),
- quadraticEqn(- x*x, c*x + a*d, d*x).roots()
+ linearEqn(-a, x).roots(),
+ quadraticEqn(-x*x, c*x + a*d, d*x).roots()
);
}
diff --git a/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.H b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.H
index f8ba53e750e920b1a2c4f1f704b94e17d67fb093..3e90378981678023713421bbd4cbeb5124372971 100644
--- a/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.H
+++ b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.H
@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2017 OpenFOAM Foundation
+ Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -27,7 +28,67 @@ Class
Foam::cubicEqn
Description
- Cubic equation of the form a*x^3 + b*x^2 + c*x + d = 0
+ Container to encapsulate various operations for
+ cubic equation of the forms with real coefficients:
+
+ \f[
+ a*x^3 + b*x^2 + c*x + d = 0
+ x^3 + B*x^2 + C*x + D = 0
+ \f]
+
+ The following two substitutions into the above forms are used:
+
+ \f[
+ x = t - B/3
+ t = w - P/3/w
+ \f]
+
+ This reduces the problem to a quadratic in \c w^3:
+
+ \f[
+ w^6 + Q*w^3 - P = 0
+ \f]
+
+ where \c Q and \c P are given in the code.
+
+ The properties of the cubic can be related to the properties of this
+ quadratic in \c w^3.
+
+ If the quadratic eqn has two identical real roots at zero, three identical
+ real roots exist in the cubic eqn.
+
+ If the quadratic eqn has two identical real roots at non-zero, two identical
+ and one distinct real roots exist in the cubic eqn.
+
+ If the quadratic eqn has two complex roots (a complex conjugate pair),
+ three distinct real roots exist in the cubic eqn.
+
+ If the quadratic eqn has two distinct real roots, one real root and two
+ complex roots (a complex conjugate pair) exist in the cubic eqn.
+
+ The quadratic eqn is solved for the most numerically accurate value
+ of \c w^3. See the \link quadraticEqn.H \endlink for details on how to
+ pick a value. This single value of \c w^3 can yield up to three cube
+ roots for \c w, which relate to the three solutions for \c x.
+
+ Only a single root, or pair of conjugate roots, is directly evaluated; the
+ one, or ones with the lowest relative numerical error. Root identities are
+ then used to recover the remaining roots, possibly utilising a quadratic
+ and/or linear solution. This seems to be a good way of maintaining the
+ accuracy of roots at very different magnitudes.
+
+ Reference:
+ \verbatim
+ Kahan's algo. to compute 'p' using fused multiply-adds (tag:JML):
+ Jeannerod, C. P., Louvet, N., & Muller, J. M. (2013).
+ Further analysis of Kahan's algorithm for the accurate
+ computation of 2× 2 determinants.
+ Mathematics of Computation, 82(284), 2245-2264.
+ DOI:10.1090/S0025-5718-2013-02679-8
+ \endverbatim
+
+See also
+ Test-cubicEqn.C
SourceFiles
cubicEqnI.H
@@ -79,28 +140,37 @@ public:
// Member Functions
- // Access
+ // Access
+
+ inline scalar a() const;
+ inline scalar b() const;
+ inline scalar c() const;
+ inline scalar d() const;
- inline scalar a() const;
- inline scalar b() const;
- inline scalar c() const;
- inline scalar d() const;
+ inline scalar& a();
+ inline scalar& b();
+ inline scalar& c();
+ inline scalar& d();
- inline scalar& a();
- inline scalar& b();
- inline scalar& c();
- inline scalar& d();
+ // Evaluate
- //- Evaluate at x
+ //- Evaluate the cubic equation at x
inline scalar value(const scalar x) const;
- //- Evaluate the derivative at x
+ //- Evaluate the derivative of the cubic equation at x
inline scalar derivative(const scalar x) const;
- //- Estimate the error of evaluation at x
+ //- Estimate the error of evaluation of the cubic equation at x
inline scalar error(const scalar x) const;
- //- Get the roots
+ //- Return the roots of the cubic equation with no particular order
+ // if discriminant > 0: return three distinct real roots
+ // if discriminant < 0: return one real root and one complex root
+ // (one member of the complex conjugate pair)
+ // if discriminant = 0: return two identical roots,
+ // and one distinct real root
+ // if identical zero Hessian: return three identical real roots
+ // where the discriminant = - 4p^3 - 27q^2.
Roots<3> roots() const;
};
diff --git a/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqn.H b/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqn.H
index adb6f39591a564a0982b0a9632d748aefa6908c6..c898aef1ab1566ed2ccfa3ee678f33f5b3a62d68 100644
--- a/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqn.H
+++ b/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqn.H
@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2017 OpenFOAM Foundation
+ Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -27,7 +28,16 @@ Class
Foam::linearEqn
Description
- Linear equation of the form a*x + b = 0
+ Container to encapsulate various operations for
+ linear equation of the forms with real coefficients:
+
+ \f[
+ a*x + b = 0
+ x + B = 0
+ \f]
+
+See also
+ Test-linearEqn.C
SourceFiles
linearEqnI.H
@@ -72,24 +82,26 @@ public:
// Member Functions
- // Access
+ // Access
+
+ inline scalar a() const;
+ inline scalar b() const;
- inline scalar a() const;
- inline scalar b() const;
+ inline scalar& a();
+ inline scalar& b();
- inline scalar& a();
- inline scalar& b();
+ // Evaluate
- //- Evaluate at x
+ //- Evaluate the linear equation at x
inline scalar value(const scalar x) const;
- //- Evaluate the derivative at x
+ //- Evaluate the derivative of the linear equation at x
inline scalar derivative(const scalar x) const;
- //- Estimate the error of evaluation at x
+ //- Estimate the error of evaluation of the linear equation at x
inline scalar error(const scalar x) const;
- //- Get the roots
+ //- Return the real root of the linear equation
inline Roots<1> roots() const;
};
diff --git a/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqnI.H b/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqnI.H
index 29b2013341f4157e83a74caa9323c8983df778fd..0f2fda08144847011443ae12f950778f8831919b 100644
--- a/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqnI.H
+++ b/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqnI.H
@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2017 OpenFOAM Foundation
+ Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -90,29 +91,19 @@ inline Foam::scalar Foam::linearEqn::error(const scalar x) const
inline Foam::Roots<1> Foam::linearEqn::roots() const
{
- /*
-
- This function solves a linear equation of the following form:
-
- a*x + b = 0
- x + B = 0
-
- */
-
const scalar a = this->a();
const scalar b = this->b();
- if (a == 0)
+ if (mag(a) < VSMALL)
{
return Roots<1>(roots::nan, 0);
}
-
- if (mag(b/VGREAT) >= mag(a))
+ else if (mag(b/VGREAT) >= mag(a))
{
return Roots<1>(sign(a) == sign(b) ? roots::negInf : roots::posInf, 0);
}
- return Roots<1>(roots::real, - b/a);
+ return Roots<1>(roots::real, -b/a);
}
diff --git a/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.C b/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.C
index ca6c3cfd3b377d8777c0b3e984d679a0f7af243c..910f9375ce6899989e141a76e9cff27d35355480 100644
--- a/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.C
+++ b/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.C
@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2017 OpenFOAM Foundation
+ Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -32,58 +33,42 @@ License
Foam::Roots<2> Foam::quadraticEqn::roots() const
{
- /*
-
- This function solves a quadraticEqn equation of the following form:
-
- a*x^2 + b*x + c = 0
- x^2 + B*x + C = 0
-
- The quadraticEqn formula is as follows:
-
- x = - B/2 +- sqrt(B*B - 4*C)/2
-
- If the sqrt generates a complex number, this provides the result. If not
- then the real root with the smallest floating point error is calculated.
-
- x0 = - B/2 - sign(B)*sqrt(B*B - 4*C)/2
-
- The other root is the obtained using an identity.
-
- x1 = C/x0
-
- */
-
const scalar a = this->a();
const scalar b = this->b();
const scalar c = this->c();
- if (a == 0)
+ // Check the leading term in the quadratic eqn exists
+ if (mag(a) < VSMALL)
{
return Roots<2>(linearEqn(b, c).roots(), roots::nan, 0);
}
- // This is assumed not to over- or under-flow. If it does, all bets are off.
- const scalar disc = b*b/4 - a*c;
+ // (JLM:p. 2246) [discriminant = b*b/4 - a*c]
+ const scalar w = a*c;
+ const scalar numDiscr = fma(-a, c, w) + fma(b, b/4, -w);
+ const scalar discr = (mag(numDiscr) > SMALL) ? numDiscr : 0;
- // How many roots of what types are available?
- const bool oneReal = disc == 0;
- const bool twoReal = disc > 0;
- //const bool twoComplex = disc < 0;
+ // Find how many roots of what types are available
+ const bool twoReal = discr > 0;
+ const bool twoComplex = discr < 0;
+ //const bool oneReal = discr == 0;
- if (oneReal)
+ if (twoReal)
{
- const Roots<1> r = linearEqn(a, b/2).roots();
- return Roots<2>(r, r);
+ // (F:Exp. 8.9)
+ const scalar x = -b/2 - sign(b)*sqrt(discr);
+ return Roots<2>(linearEqn(-a, x).roots(), linearEqn(-x, c).roots());
}
- else if (twoReal)
+ else if (twoComplex)
{
- const scalar x = - b/2 - sign(b)*sqrt(disc);
- return Roots<2>(linearEqn(- a, x).roots(), linearEqn(- x, c).roots());
+ const Roots<1> xRe(roots::type::complex, -b/2/a);
+ const Roots<1> xIm(roots::type::complex, sign(b)*sqrt(mag(discr))/a);
+ return Roots<2>(xRe, xIm);
}
- else // if (twoComplex)
+ else // (oneReal)
{
- return Roots<2>(roots::complex, 0);
+ const Roots<1> r(linearEqn(a, b/2).roots());
+ return Roots<2>(r, r);
}
}
diff --git a/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.H b/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.H
index 8e6d6dd9fc0ef6360dbf01a7953432985b614832..c68da9b9ec4d9014699e535f0692b02a95092a10 100644
--- a/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.H
+++ b/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.H
@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2017 OpenFOAM Foundation
+ Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -27,7 +28,43 @@ Class
Foam::quadraticEqn
Description
- Quadratic equation of the form a*x^2 + b*x + c = 0
+ Container to encapsulate various operations for
+ quadratic equation of the forms with real coefficients:
+
+ \f[
+ a*x^2 + b*x + c = 0
+ x^2 + B*x + C = 0
+ \f]
+
+ The expressions for the roots of \c quadraticEqn are as follows:
+
+ \f[
+ x1 = - (b + sign(b) sqrt(b^2 - 4ac)/(2*a))
+
+ x2 = c/(a*x1)
+ \f]
+
+ where \c (b^2 - 4ac) is evaluated by fused multiply-adds to avoid
+ detrimental cancellation.
+
+ Reference:
+ \verbatim
+ Cancellation-avoiding quadratic formula (tag:F):
+ Ford, W. (2014).
+ Numerical linear algebra with applications: Using MATLAB.
+ London: Elsevier/Academic Press.
+ DOI:10.1016/C2011-0-07533-6
+
+ Kahan's algo. to compute 'b^2-a*c' using fused multiply-adds (tag:JML):
+ Jeannerod, C. P., Louvet, N., & Muller, J. M. (2013).
+ Further analysis of Kahan's algorithm for the accurate
+ computation of 2× 2 determinants.
+ Mathematics of Computation, 82(284), 2245-2264.
+ DOI:10.1090/S0025-5718-2013-02679-8
+ \endverbatim
+
+See also
+ Test-quadraticEqn.C
SourceFiles
quadraticEqnI.H
@@ -73,26 +110,31 @@ public:
// Member Functions
- // Access
+ // Access
+
+ inline scalar a() const;
+ inline scalar b() const;
+ inline scalar c() const;
- inline scalar a() const;
- inline scalar b() const;
- inline scalar c() const;
+ inline scalar& a();
+ inline scalar& b();
+ inline scalar& c();
- inline scalar& a();
- inline scalar& b();
- inline scalar& c();
+ // Evaluate
- //- Evaluate at x
+ //- Evaluate the quadratic equation at x
inline scalar value(const scalar x) const;
- //- Evaluate the derivative at x
+ //- Evaluate the derivative of the quadratic equation at x
inline scalar derivative(const scalar x) const;
- //- Estimate the error of evaluation at x
+ //- Estimate the error of evaluation of the quadratic equation at x
inline scalar error(const scalar x) const;
- //- Get the roots
+ //- Return the roots of the quadratic equation with no particular order
+ // if discriminant > 0: return two distinct real roots
+ // if discriminant < 0: return one of the complex conjugate-pair roots
+ // otherwise : return two identical real roots
Roots<2> roots() const;
};