diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj
index 20d19ca8..7f874aff 100644
--- a/src/Numerics/Numerics.csproj
+++ b/src/Numerics/Numerics.csproj
@@ -191,6 +191,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 ebe5f2ce..1bd33bee 100644
--- a/src/UnitTests/UnitTests.csproj
+++ b/src/UnitTests/UnitTests.csproj
@@ -378,6 +378,7 @@
+