Roots of polynomials with Complex Coefficients

src/Numerics.Tests/RootFindingTests/FindRootsTest.cs

@@ -27,9 +27,10 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
+using NUnit.Framework;
 using System;
+using System.Linq;
 using Complex = System.Numerics.Complex;
-using NUnit.Framework;
 
 namespace MathNet.Numerics.Tests.RootFindingTests
 {
@@ -210,5 +211,44 @@ public void QuadraticExpected(double c, double b, double a, double x1R, double x
             AssertComplexEqual(Complex.Zero, c + b*x1 + a*x1*x1, 1e-14);
             AssertComplexEqual(Complex.Zero, c + b*x2 + a*x2*x2, 1e-14);
         }
+
+        [Test]
+        public void RootsTest()
+        {           
+            var c1 = new Complex[] { new Complex(1, 1), new Complex (1,1), new Complex (1,1) };
+
+            var r1 = FindRoots.Roots( new Complex[] { c1[0] });
+            Assert.AreEqual(Array.Empty<Complex>(), r1);
+
+            var r2 = FindRoots.Roots(new Complex[] { c1[0], c1[1] });
+            Assert.AreEqual (Array.Empty<Complex>(), r2);
+
+            //expected imaginary component = Math.Sqrt(3) / 2
+            var expected_1 = new Complex[] { new Complex(-0.5, -0.8660254037844386), new Complex(-0.5, 0.8660254037844386) };
+            TestRootsEqual(c1, expected_1);
+
+            //Expected root values calculated with cube root forumla using System.Math functions
+            var c2 = new Complex[] { new Complex(-2, 3), new Complex(1, 0), new Complex(0, 1), new Complex(2, 1)};
+            var expected_2 = new Complex[] { new Complex(0.828084722641885, -0.6508989159005254), new Complex(-0.9759670820130891, -0.900661591057162), new Complex(-0.052117640628795536, 1.1515605069576873) };
+            TestRootsEqual (c2, expected_2);
+        }
+
+        static void TestRootsEqual(Complex[] x, Complex[] eIn)
+        {
+            var tol = 1e-10;
+            var r0 = FindRoots.Roots(x);
+
+            var e = eIn.OrderBy(v => v.Imaginary).ToArray();
+            var r = r0.OrderBy(v => v.Imaginary).ToArray();
+
+            Assert.IsNotEmpty(r);
+            Assert.AreEqual(e.Length, r.Length, "Length mismatch");
+            for (int k = 0; k < r.Length; k++)
+            {
+                var msg = String.Format("At k={0}", k);
+                Assert.AreEqual(e[k].Real, r[k].Real, tol, msg);
+                Assert.AreEqual(e[k].Imaginary, r[k].Imaginary, tol, msg);
+            }           
+        }
     }
 }

src/Numerics/FindRoots.cs

@@ -27,8 +27,12 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
-using System;
+using MathNet.Numerics.LinearAlgebra;
+using MathNet.Numerics.LinearAlgebra.Complex;
+using MathNet.Numerics.LinearAlgebra.Factorization;
 using MathNet.Numerics.RootFinding;
+using System;
+using System.Linq;
 using Complex = System.Numerics.Complex;
 
 namespace MathNet.Numerics
@@ -118,6 +122,39 @@ public static Complex[] Polynomial(double[] coefficients)
             return new Polynomial(coefficients).Roots();
         }
 
+        /// <summary>
+        /// Find all roots of a polynomial with complex coefficients by calculating the characteristic polynomial of the companion matrix
+        /// </summary>
+        /// <param name="coefficients">The coefficients of the polynomial in ascending order, e.g. new Complex[] {new (5, 0), new (0, 2), new (2,3)} = "5 + 2i x^1 + (2 + 3i) x^2"</param>
+        /// <returns>The roots of the polynomial</returns>
+
+        public static Complex[] Roots(Complex[] coefficients)
+        {
+            int n = coefficients.Length - 1;
+            if (n < 2)
+            {
+                return Array.Empty<Complex>();
+            }
+
+            // Negate, and normalize (scale such that the polynomial becomes monic)
+            Complex aN = coefficients[n];
+            Complex[] p = new Complex[n];
+            for (int i = n - 1; i >= 0; i--)
+            {
+                p[i] = -coefficients[i] / aN;
+            }
+
+            DenseMatrix A0 = DenseMatrix.CreateDiagonal(n - 1, n - 1, 1.0);
+            DenseMatrix A = new DenseMatrix(n);
+
+            A.SetSubMatrix(1, 0, A0);
+            A.SetRow(0, p.Reverse().ToArray());
+
+            Evd<Complex> eigen = A.Evd(Symmetricity.Asymmetric);
+
+            return eigen.EigenValues.ToArray();
+        }
+
         /// <summary>
         /// Find all roots of a polynomial by calculating the characteristic polynomial of the companion matrix
         /// </summary>