From 09c89cc4005bef3dd1c3637b504954ebc3cb3fe7 Mon Sep 17 00:00:00 2001 From: Candy Chiu Date: Mon, 20 Jan 2014 20:04:41 -0500 Subject: [PATCH] Added class Cubic to find the real roots of a cubic polynomial. Added unit test class CubicTest. --- src/Numerics/Numerics.csproj | 1 + src/Numerics/RootFinding/Cubic.cs | 96 +++++++++++++++++++++ src/UnitTests/RootFindingTests/CubicTest.cs | 61 +++++++++++++ src/UnitTests/UnitTests.csproj | 1 + 4 files changed, 159 insertions(+) create mode 100644 src/Numerics/RootFinding/Cubic.cs create mode 100644 src/UnitTests/RootFindingTests/CubicTest.cs diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 31785ea8..b3ac67ee 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -189,6 +189,7 @@ + diff --git a/src/Numerics/RootFinding/Cubic.cs b/src/Numerics/RootFinding/Cubic.cs new file mode 100644 index 00000000..9d21ae7d --- /dev/null +++ b/src/Numerics/RootFinding/Cubic.cs @@ -0,0 +1,96 @@ +using System; + +#if !NOSYSNUMERICS + using Complex = System.Numerics.Complex; +#endif + +namespace MathNet.Numerics.RootFinding +{ + /// + /// Finds roots to the cubic equation x^3 + a2*x^2 + a1*x + a0 = 0 + /// Implements the cubic formula in http://mathworld.wolfram.com/CubicFormula.html + /// + public static class Cubic + { + // D = Q^3 + R^2 is the polynomial discriminant. + // D > 0, 1 real root + // D = 0, 3 real roots, at least two are equal + // D < 0, 3 real and unequal roots + + /// + /// Q and R are transformed variables. + /// + private static void QR(double a2, double a1, double a0, ref double Q, ref double R) + { + Q = (3 * a1 - a2 * a2)/9.0; + R = (9.0 * a2 * a1 - 27 * a0 - 2 * a2 * a2 * a2)/54.0; + } + + /// + /// n^(1/3) - work around a negative double raised to (1/3) + /// + private static double PowThird(double n) + { + return Math.Pow(Math.Abs(n), 1d / 3d) * Math.Sign(n); + } + + public static Tuple RealRoots(double a2, double a1, double a0) + { + var Q = double.NaN; + var R = double.NaN; + QR(a2, a1, a0, ref Q, ref R); + + var Q3 = Q * Q * Q; + var D = Q3 + R * R; + var shift = -a2 / 3d; + + double x1 = double.NaN; + double x2 = double.NaN; + double x3 = double.NaN; + + // when D >= 0, use eqn (54)-(56) where S and T are real + if (D >= 0) + { + double sqrtD = Math.Pow(D, 0.5); + double S = PowThird(R + sqrtD); + double T = PowThird(R - sqrtD); + x1 = shift + (S + T); + if (D == 0) + x2 = shift - S; + } + // 3 real roots, use eqn (70)-(73) to calculate the real roots + else + { + double theta = Math.Acos(R / Math.Sqrt(-Q3)); + x1 = 2d * Math.Sqrt(-Q) * Math.Cos(theta / 3.0) + shift; + x2 = 2d * Math.Sqrt(-Q) * Math.Cos((theta + 2.0 * Constants.Pi) / 3d) + shift; + x3 = 2d * Math.Sqrt(-Q) * Math.Cos((theta - 2.0 * Constants.Pi) / 3d) + shift; + } + return Tuple.Create(x1, x2, x3); + } + + public static Tuple Roots(double a2, double a1, double a0) + { + // use eqn (54)-(56) + var Q = double.NaN; + var R = double.NaN; + QR(a2, a1, a0, ref Q, ref R); + + var D = Q * Q * Q + R * R; + + var rootD = Complex.Sqrt(D); + var S = Complex.Pow(R + rootD, 1d / 3d); + var T = Complex.Pow(R - rootD, 1d / 3d); + var shift = -a2 / 3d; + var sharedI = 0.5 * Complex.ImaginaryOne * Math.Sqrt(3) * (S - T); + + var x1 = shift + (S + T); + var x2 = shift - 0.5 * (S + T); + var x3 = x2; + x2 += sharedI; + x3 -= sharedI; + + return Tuple.Create(x1, x2, x3); + } + } +} diff --git a/src/UnitTests/RootFindingTests/CubicTest.cs b/src/UnitTests/RootFindingTests/CubicTest.cs new file mode 100644 index 00000000..c8d71a3f --- /dev/null +++ b/src/UnitTests/RootFindingTests/CubicTest.cs @@ -0,0 +1,61 @@ +using System; +using MathNet.Numerics.RootFinding; +using NUnit.Framework; + +namespace MathNet.Numerics.UnitTests.RootFindingTests +{ + [TestFixture, Category("RootFinding")] + internal class CubicTest + { + [Test] + public void RealRoots_SingleReal() + { + var a0 = 1d; + var a1 = 5d; + var a2 = 2d; + + var roots = Cubic.RealRoots(a2, a1, a0); + var root = roots.Item1; + var funcValue = root * (root * (root + a2) + a1) + a0; + Assert.That(funcValue, Is.EqualTo(0).Within(1e-14)); + Assert.AreEqual(double.NaN, roots.Item2); + Assert.AreEqual(double.NaN, roots.Item3); + } + + [Test] + public void RealRoots_DoubleReal1() + { + var a0 = -2d; + var a1 = 5d; + var a2 = -4d; + var roots = Cubic.RealRoots(a2, a1, a0); + Assert.That(roots.Item1, Is.EqualTo(2).Within(1e-14)); + Assert.That(roots.Item2, Is.EqualTo(1).Within(1e-14)); + Assert.AreEqual(double.NaN, roots.Item3); + } + + [Test] + public void RealRoots_DoubleReal2() + { + var a0 = -4d; + var a1 = 8d; + var a2 = -5d; + var roots = Cubic.RealRoots(a2, a1, a0); + Assert.That(roots.Item1, Is.EqualTo(1).Within(1e-14)); + Assert.That(roots.Item2, Is.EqualTo(2).Within(1e-14)); + Assert.AreEqual(double.NaN, roots.Item3); + } + + [Test] + public void RealRoots_TripleReal() + { + var a0 = 6d; + var a1 = -5d; + var a2 = -2d; + var roots = Cubic.RealRoots(a2, a1, a0); + Assert.That(roots.Item1, Is.EqualTo(3).Within(1e-14)); + Assert.That(roots.Item2, Is.EqualTo(-2).Within(1e-14)); + Assert.That(roots.Item3, Is.EqualTo(1).Within(1e-14)); + } + } +} diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 54104d8d..b54c9a3f 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -376,6 +376,7 @@ +