From 9a641040f9119dde45f7f8f5adc1d02a580d79a1 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sun, 14 Oct 2018 20:51:38 +0200 Subject: [PATCH] Polynomial: reimplement arithmetics in immutable way --- src/Numerics.Tests/PolynomialTests.cs | 57 +++---- src/Numerics/Polynomial.cs | 206 +++++++++++++------------- 2 files changed, 128 insertions(+), 135 deletions(-) diff --git a/src/Numerics.Tests/PolynomialTests.cs b/src/Numerics.Tests/PolynomialTests.cs index 2a45d7f0..0d559c15 100644 --- a/src/Numerics.Tests/PolynomialTests.cs +++ b/src/Numerics.Tests/PolynomialTests.cs @@ -66,8 +66,8 @@ namespace MathNet.Numerics.UnitTests } else { - Assert.AreEqual(expected.Length, p_res.CoefficientCount, "length mismatch"); - for (int k = 0; k < p_res.CoefficientCount; k++) + Assert.AreEqual(expected.Length, p_res.Coefficients.Length, "length mismatch"); + for (int k = 0; k < p_res.Coefficients.Length; k++) { Assert.AreEqual(expected[k], p_res.Coefficients[k], "idx: " + k + " mismatch"); } @@ -89,8 +89,8 @@ namespace MathNet.Numerics.UnitTests } else { - Assert.AreEqual(expected.Length, p_res.CoefficientCount, "length mismatch"); - for (int k = 0; k < p_res.CoefficientCount; k++) + Assert.AreEqual(expected.Length, p_res.Coefficients.Length, "length mismatch"); + for (int k = 0; k < p_res.Coefficients.Length; k++) { Assert.AreEqual(expected[k], p_res.Coefficients[k], "idx: " + k + " mismatch"); } @@ -125,8 +125,8 @@ namespace MathNet.Numerics.UnitTests p_res.Trim(); p_tar.Trim(); - Assert.AreEqual(p_tar.CoefficientCount, p_res.CoefficientCount, "length mismatch"); - for (int k = 0; k < p_res.CoefficientCount; k++) + Assert.AreEqual(p_tar.Coefficients.Length, p_res.Coefficients.Length, "length mismatch"); + for (int k = 0; k < p_res.Coefficients.Length; k++) { Assert.AreEqual(p_tar.Coefficients[k], p_res.Coefficients[k], msg); } @@ -135,7 +135,7 @@ namespace MathNet.Numerics.UnitTests } [Test] - public void SubstractTest() + public void SubtractTest() { for (int i = 0; i < 5; i++) { @@ -156,14 +156,14 @@ namespace MathNet.Numerics.UnitTests var p1 = new Polynomial(c1); var p2 = new Polynomial(c2); - var p_res = Polynomial.Substract(p1, p2); + var p_res = Polynomial.Subtract(p1, p2); var p_tar = new Polynomial(tgt); p_res.Trim(); p_tar.Trim(); - Assert.AreEqual(p_tar.CoefficientCount, p_res.CoefficientCount, "length mismatch"); - for (int k = 0; k < p_res.CoefficientCount; k++) + Assert.AreEqual(p_tar.Coefficients.Length, p_res.Coefficients.Length, "length mismatch"); + for (int k = 0; k < p_res.Coefficients.Length; k++) { Assert.AreEqual(p_tar.Coefficients[k], p_res.Coefficients[k], msg); } @@ -198,8 +198,8 @@ namespace MathNet.Numerics.UnitTests p_res.Trim(); p_tar.Trim(); - Assert.AreEqual(p_tar.CoefficientCount, p_res.CoefficientCount, "length mismatch"); - for (int k = 0; k < p_res.CoefficientCount; k++) + Assert.AreEqual(p_tar.Coefficients.Length, p_res.Coefficients.Length, "length mismatch"); + for (int k = 0; k < p_res.Coefficients.Length; k++) { Assert.AreEqual(p_tar.Coefficients[k], p_res.Coefficients[k], msg); } @@ -256,14 +256,14 @@ namespace MathNet.Numerics.UnitTests var p11 = new Polynomial(2.0d); var p21 = new Polynomial(2.0d); var tpl1 = Polynomial.DivideLong(p11, p21); - testEqual(new double[] { 1.0 }, tpl1.Item1); - testEqual(new double[] { 0.0 }, tpl1.Item2); + TestEqual(new double[] { 1.0 }, tpl1.Item1); + TestEqual(new double[] { 0.0 }, tpl1.Item2); var p12 = new Polynomial(new double[] { 2.0d, 2.0d }); var p22 = new Polynomial(2.0d); var tpl2 = Polynomial.DivideLong(p12, p22); - testEqual(new double[] { 1.0, 1.0 }, tpl2.Item1); - testEqual(new double[] { 0.0 }, tpl2.Item2); + TestEqual(new double[] { 1.0, 1.0 }, tpl2.Item1); + TestEqual(new double[] { 0.0 }, tpl2.Item2); for (int i = 0; i < 5; i++) { @@ -286,7 +286,7 @@ namespace MathNet.Numerics.UnitTests var pres = (pquo * pi) + prem; pres.Trim(); - testEqual(pres, tgt, msg); + TestEqual(pres, tgt, msg); } } } @@ -314,23 +314,23 @@ namespace MathNet.Numerics.UnitTests var x_2 = new double[] { -1.0, 1.0 }; var expected_2 = new List(); expected_2.Add(new Complex(1.0, 0.0)); - testEqual(x_2, expected_2); + TestEqual(x_2, expected_2); var x_3 = new double[] { -1.0, 0.0, 1.0 }; var expected_3 = new List(); expected_3.Add(new Complex(1.0, 0.0)); expected_3.Add(new Complex(-1.0, 0.0)); - testEqual(x_3, expected_3); + TestEqual(x_3, expected_3); var x_4 = new double[] { -1.0, -0.33333333333333337, 0.33333333333333326, 1.0 }; var expected_4 = new List(); expected_4.Add(new Complex(0.9999999999999996, 0.0)); expected_4.Add(new Complex(-0.6666666666666666, 0.7453559924999296)); expected_4.Add(new Complex(-0.6666666666666666, -0.7453559924999296)); - testEqual(x_4, expected_4); + TestEqual(x_4, expected_4); } - private void testEqual(double[] x, List eIn) + static void TestEqual(double[] x, List eIn) { var tol = 1e-10; var r0 = new Polynomial(x).GetRoots().ToList(); @@ -348,7 +348,7 @@ namespace MathNet.Numerics.UnitTests } } - private void testEqual(double[] p_tar, double[] p_res, string msg = null) + static void TestEqual(double[] p_tar, double[] p_res, string msg = null) { Assert.AreEqual(p_tar.Length, p_res.Length, "length mismatch"); for (int k = 0; k < p_res.Length; k++) @@ -357,19 +357,20 @@ namespace MathNet.Numerics.UnitTests } } - private void testEqual(double[] p_tar, Polynomial p_res, string msg = null) + static void TestEqual(double[] p_tar, Polynomial p_res, string msg = null) { - Assert.AreEqual(p_tar.Length, p_res.CoefficientCount, "length mismatch"); - for (int k = 0; k < p_res.CoefficientCount; k++) + Assert.AreEqual(p_tar.Length, p_res.Coefficients.Length, "length mismatch"); + for (int k = 0; k < p_res.Coefficients.Length; k++) { Assert.AreEqual(p_tar[k], p_res.Coefficients[k], msg); } } - private void testEqual(Polynomial p_tar, Polynomial p_res, string msg = null) + + static void TestEqual(Polynomial p_tar, Polynomial p_res, string msg = null) { - Assert.AreEqual(p_tar.CoefficientCount, p_res.CoefficientCount, "length mismatch"); - for (int k = 0; k < p_res.CoefficientCount; k++) + Assert.AreEqual(p_tar.Coefficients.Length, p_res.Coefficients.Length, "length mismatch"); + for (int k = 0; k < p_res.Coefficients.Length; k++) { Assert.AreEqual(p_tar.Coefficients[k], p_res.Coefficients[k], msg); } diff --git a/src/Numerics/Polynomial.cs b/src/Numerics/Polynomial.cs index 1ed00007..d79056fd 100644 --- a/src/Numerics/Polynomial.cs +++ b/src/Numerics/Polynomial.cs @@ -17,24 +17,13 @@ namespace MathNet.Numerics /// /// The coefficients of the polynomial in a /// - public double[] Coefficients { get; set; } + public double[] Coefficients { get; private set; } /// /// Only needed for the ToString method /// public string VarName = "x^"; - /// - /// Length of Polynomial (max element + 1) e.G x^5 highest element, will give Length = 6 - /// - public int CoefficientCount - { - get - { - return (Coefficients == null ? 0 : Coefficients.Length); - } - } - /// /// Degree of the polynomial, i.e. the largest monomial exponent. For example, the degree of x^2+x^5 is 5. /// The null-polynomial returns degree -1 because the correct degree, negative infinity, cannot be represented by integers. @@ -115,17 +104,25 @@ namespace MathNet.Numerics /// public void Trim() { - if (CoefficientCount == 1) + if (Coefficients.Length == 1) + { return; + } - int i = CoefficientCount - 1; + int i = Coefficients.Length - 1; while (i >= 0 && Coefficients[i] == 0.0) + { i--; + } if (i < 0) + { Coefficients = new[] { 0.0 }; + } else if (i == 0) + { Coefficients = new[] { Coefficients[0] }; + } else { var hold = new double[i+1]; @@ -304,7 +301,7 @@ namespace MathNet.Numerics /// resulting Polynomial public static Polynomial operator -(Polynomial a, Polynomial b) { - return Substract(a, b); + return Subtract(a, b); } /// @@ -368,128 +365,147 @@ namespace MathNet.Numerics } /// - /// pointwise division of two Polynomials + /// Point-wise division of two Polynomials /// - /// left Polynomial - /// right Polynomial - /// resulting Polynomial - public static Polynomial DividePointwise(Polynomial a, Polynomial b) + /// Left Polynomial + /// Right Polynomial + /// Resulting Polynomial + public static Polynomial PointwiseDivide(Polynomial a, Polynomial b) { - var aa = a.Clone() as Polynomial; - var bb = b.Clone() as Polynomial; + var ac = a.Coefficients; + var bc = b.Coefficients; - if (aa.Coefficients.Length != bb.Coefficients.Length) + var degree = a.Degree; + var result = new double[degree + 1]; + + var commonLength = Math.Min(Math.Min(ac.Length, bc.Length), result.Length); + for (int i = 0; i < commonLength; i++) { - MakeSameLength(ref aa, ref bb); + result[i] = ac[i] / bc[i]; } - int n = aa.Coefficients.Length; - double[] res = new double[aa.Coefficients.Length]; - for (int ii = 0; ii < n; ii++) + for (int i = commonLength; i < result.Length; i++) { - res[ii] = aa.Coefficients[ii] / bb.Coefficients[ii]; + result[i] = ac[i] / 0.0; } - return new Polynomial(res); + return new Polynomial(result); } /// - /// pointwise multiplication of two Polynomials + /// Point-wise multiplication of two Polynomials /// - /// left Polynomial - /// right Polynomial - /// resulting Polynomial - public static Polynomial MultiplyPointwise(Polynomial a, Polynomial b) + /// Left Polynomial + /// Right Polynomial + /// Resulting Polynomial + public static Polynomial PointwiseMultiply(Polynomial a, Polynomial b) { - var aa = a.Clone() as Polynomial; - var bb = b.Clone() as Polynomial; - - if (aa.Coefficients.Length != bb.Coefficients.Length) - { - MakeSameLength(ref aa, ref bb); - } + var ac = a.Coefficients; + var bc = b.Coefficients; - int n = aa.Coefficients.Length; - double[] res = new double[aa.Coefficients.Length]; - for (int ii = 0; ii < n; ii++) + var degree = Math.Min(a.Degree, b.Degree); + var result = new double[degree + 1]; + for (int i = 0; i < result.Length; i++) { - res[ii] = aa.Coefficients[ii] * bb.Coefficients[ii]; + result[i] = ac[i] * bc[i]; } - return new Polynomial(res); + return new Polynomial(result); } /// - /// Addition of two Polynomials (piecewise) + /// Addition of two Polynomials (point-wise). /// - /// left Polynomial - /// right Polynomial - /// resulting Polynomial + /// Left Polynomial + /// Right Polynomial + /// Resulting Polynomial public static Polynomial Add(Polynomial a, Polynomial b) { - var aa = a.Clone() as Polynomial; - var bb = b.Clone() as Polynomial; + var ac = a.Coefficients; + var bc = b.Coefficients; - if (aa.CoefficientCount != bb.CoefficientCount) + var degree = Math.Max(a.Degree, b.Degree); + var result = new double[degree + 1]; + + var commonLength = Math.Min(Math.Min(ac.Length, bc.Length), result.Length); + for (int i = 0; i < commonLength; i++) { - MakeSameLength(ref aa, ref bb); + result[i] = ac[i] + bc[i]; } - int n = aa.CoefficientCount; - double[] res = new double[n]; - for (int ii = 0; ii < n; ii++) + int acLength = Math.Min(ac.Length, result.Length); + for (int i = commonLength; i < acLength; i++) { - res[ii] = aa.Coefficients[ii] + bb.Coefficients[ii]; + // no need to add since only one of both applies + result[i] = ac[i]; } - return new Polynomial(res); + int bcLength = Math.Min(bc.Length, result.Length); + for (int i = commonLength; i < bcLength; i++) + { + // no need to add since only one of both applies + result[i] = bc[i]; + } + + return new Polynomial(result); } /// - /// substraction of two Polynomials (piecewise) + /// Subtraction of two Polynomials (point-wise). /// - /// left Polynomial - /// right Polynomial - /// resulting Polynomial - public static Polynomial Substract(Polynomial a, Polynomial b) + /// Left Polynomial + /// Right Polynomial + /// Resulting Polynomial + public static Polynomial Subtract(Polynomial a, Polynomial b) { - var aa = a.Clone() as Polynomial; - var bb = b.Clone() as Polynomial; + var ac = a.Coefficients; + var bc = b.Coefficients; + + var degree = Math.Max(a.Degree, b.Degree); + var result = new double[degree + 1]; + + var commonLength = Math.Min(Math.Min(ac.Length, bc.Length), result.Length); + for (int i = 0; i < commonLength; i++) + { + result[i] = ac[i] - bc[i]; + } - if (aa.CoefficientCount != bb.CoefficientCount) + int acLength = Math.Min(ac.Length, result.Length); + for (int i = commonLength; i < acLength; i++) { - MakeSameLength(ref aa, ref bb); + // no need to add since only one of both applies + result[i] = ac[i]; } - int n = aa.CoefficientCount; - double[] res = new double[n]; - for (int ii = 0; ii < n; ii++) + int bcLength = Math.Min(bc.Length, result.Length); + for (int i = commonLength; i < bcLength; i++) { - res[ii] = aa.Coefficients[ii] - bb.Coefficients[ii]; + // no need to add since only one of both applies + result[i] = -bc[i]; } - return new Polynomial(res); + return new Polynomial(result); } /// /// Division of two polynomials returning the quotient-with-remainder of the two polynomials given /// - /// left polynomial - /// right polynomial - /// a tuple holding quotient in first and remainder in second + /// Left polynomial + /// Right polynomial + /// A tuple holding quotient in first and remainder in second public static Tuple DivideLong(Polynomial a, Polynomial b) { if (a == null) - throw new ArgumentNullException("a"); + throw new ArgumentNullException(nameof(a)); if (b == null) - throw new ArgumentNullException("b"); + throw new ArgumentNullException(nameof(b)); - if (a.CoefficientCount <= 0) + if (a.Coefficients.Length <= 0) throw new ArgumentOutOfRangeException("a Degree must be greater than zero"); - if (b.CoefficientCount <= 0) + if (b.Coefficients.Length <= 0) throw new ArgumentOutOfRangeException("b Degree must be greater than zero"); - if (b.Coefficients[b.CoefficientCount-1] == 0) + if (b.Coefficients[b.Coefficients.Length - 1] == 0) throw new DivideByZeroException("b polynomial ends with zero"); var c1 = a.Coefficients.ToArray(); @@ -572,8 +588,8 @@ namespace MathNet.Numerics /// /// Division of two polynomials returning the quotient-with-remainder of the two polynomials given /// - /// right polynomial - /// a tuple holding quotient in first and remainder in second + /// Right polynomial + /// A tuple holding quotient in first and remainder in second public Tuple DivideLong(Polynomial b) { return DivideLong(this, b); @@ -643,30 +659,6 @@ namespace MathNet.Numerics return Coefficients.ToArray(); } - static void MakeSameLength(ref Polynomial a, ref Polynomial b) - { - double[] aHold = new double[a.Coefficients.Length]; - double[] bHold = new double[b.Coefficients.Length]; - Array.Copy(a.Coefficients, aHold, a.Coefficients.Length); - Array.Copy(b.Coefficients, bHold, b.Coefficients.Length); - - if (a.Coefficients.Length < b.Coefficients.Length) - { - a.Coefficients = new double[b.Coefficients.Length]; - b.Coefficients = new double[b.Coefficients.Length]; - Array.Copy(aHold, a.Coefficients, aHold.Length); - Array.Copy(bHold, b.Coefficients, bHold.Length); - } - else - { - a.Coefficients = new double[a.Coefficients.Length]; - b.Coefficients = new double[a.Coefficients.Length]; - Array.Copy(aHold, a.Coefficients, aHold.Length); - Array.Copy(bHold, b.Coefficients, bHold.Length); - } - - } - /// /// Full convolution of two arrays ///