diff --git a/src/Numerics.Tests/PolynomialTests.cs b/src/Numerics.Tests/PolynomialTests.cs index 7267d403..41721553 100644 --- a/src/Numerics.Tests/PolynomialTests.cs +++ b/src/Numerics.Tests/PolynomialTests.cs @@ -29,6 +29,7 @@ using System; using System.Collections.Generic; +using System.Diagnostics; using System.Linq; using System.Numerics; using MathNet.Numerics; @@ -123,8 +124,8 @@ namespace MathNet.Numerics.UnitTests var p_res = Polynomial.Add(p1, p2); var p_tar = new Polynomial(tgt); - p_res.CutTrailZeros(); - p_tar.CutTrailZeros(); + p_res.TrimTrailingZeros(); + p_tar.TrimTrailingZeros(); Assert.AreEqual(p_tar.Degree, p_res.Degree, "length mismatch"); for (int k = 0; k < p_res.Degree; k++) @@ -160,8 +161,8 @@ namespace MathNet.Numerics.UnitTests var p_res = Polynomial.Substract(p1, p2); var p_tar = new Polynomial(tgt); - p_res.CutTrailZeros(); - p_tar.CutTrailZeros(); + p_res.TrimTrailingZeros(); + p_tar.TrimTrailingZeros(); Assert.AreEqual(p_tar.Degree, p_res.Degree, "length mismatch"); for (int k = 0; k < p_res.Degree; k++) @@ -196,8 +197,8 @@ namespace MathNet.Numerics.UnitTests var p_res = p1 * p2; var p_tar = new Polynomial(tgt); - p_res.CutTrailZeros(); - p_tar.CutTrailZeros(); + p_res.TrimTrailingZeros(); + p_tar.TrimTrailingZeros(); Assert.AreEqual(p_tar.Degree, p_res.Degree, "length mismatch"); for (int k = 0; k < p_res.Degree; k++) @@ -208,6 +209,82 @@ namespace MathNet.Numerics.UnitTests } } + + public void DivideLongTestScalar(Tuple inVals, Tuple expectedVals) + { + var p1 = new Polynomial(1.0d); + var p2 = new Polynomial(new double[0]); + var tpl = Polynomial.DivideLong(p1, p2); + + + } + + [Test] + public void DivideLongTest() + { + Assert.Throws(typeof(ArgumentOutOfRangeException), () => + { + var p1 = new Polynomial(1.0d); + var p2 = new Polynomial(new double[0]); + var tpl = Polynomial.DivideLong(p1, p2); + }); + Assert.Throws(typeof(ArgumentOutOfRangeException), () => + { + var p1 = new Polynomial(1.0d); + var p2 = new Polynomial(new double[0]); + var tpl = Polynomial.DivideLong(p2, p1); + }); + Assert.Throws(typeof(ArgumentOutOfRangeException), () => + { + var p1 = new Polynomial(new double[0]); + var p2 = new Polynomial(new double[0]); + var tpl = Polynomial.DivideLong(p2, p1); + }); + Assert.Throws(typeof(DivideByZeroException), () => + { + var p1 = new Polynomial(1.0d); + var p2 = new Polynomial(0.0d); + var tpl = Polynomial.DivideLong(p1, p2); + }); + + 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); + + 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); + + for (int i = 0; i < 5; i++) + { + for (int j = 0; j < 5; j++) + { + var msg = String.Format("At i={0}, j={1}", i, j); + var ci = new double[Math.Max(2, i+1)]; + var cj = new double[Math.Max(2, j+1)]; + ci[ci.Length - 1] = 2; + ci[ci.Length - 2] = 1; + cj[cj.Length - 1] = 2; + cj[cj.Length - 2] = 1; + + var pi = new Polynomial(ci); + var pj = new Polynomial(cj); + var tgt = Polynomial.Add(pi, pj); + var tpl3 = Polynomial.DivideLong(tgt, pi); + var pquo = tpl3.Item1; + var prem = tpl3.Item2; + var pres = (pquo * pi) + prem; + testEqual(pres, tgt, msg); + } + } + } + + + [Test] public void GetRootsTest() { @@ -263,7 +340,33 @@ namespace MathNet.Numerics.UnitTests Assert.AreEqual(e[k].Real, r[k].Real, tol, msg); Assert.AreEqual(e[k].Imaginary, r[k].Imaginary, tol, msg); } + } + private 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++) + { + Assert.AreEqual(p_tar[k], p_res[k], msg); + } + } + + private void testEqual(double[] p_tar, Polynomial p_res, string msg = null) + { + Assert.AreEqual(p_tar.Length, p_res.Degree, "length mismatch"); + for (int k = 0; k < p_res.Degree; k++) + { + Assert.AreEqual(p_tar[k], p_res.Coeffs[k], msg); + } + + } + private void testEqual(Polynomial p_tar, Polynomial p_res, string msg = null) + { + Assert.AreEqual(p_tar.Degree, p_res.Degree, "length mismatch"); + for (int k = 0; k < p_res.Degree; k++) + { + Assert.AreEqual(p_tar.Coeffs[k], p_res.Coeffs[k], msg); + } } } } diff --git a/src/Numerics/Polynomial.cs b/src/Numerics/Polynomial.cs index 1da5cf91..1abd363f 100644 --- a/src/Numerics/Polynomial.cs +++ b/src/Numerics/Polynomial.cs @@ -15,11 +15,13 @@ using MathNet.Numerics.LinearAlgebra.Factorization; namespace MathNet.Numerics { /// - /// a class handlin REAL VALUED Polynomials, complex coefficients can not be handled (yet) + /// a class handling REAL VALUED Polynomials, complex coefficients can not be handled (yet) /// public class Polynomial { - + /// + /// The coefficients of the polynomial in a + /// public double[] Coeffs { get; set; } /// @@ -116,7 +118,7 @@ namespace MathNet.Numerics /// /// remove all trailing zeros, e.G before: "0.00 x^2 + 1.0 x^1 + 1.00" after: "1.0 x^1 + 1.00" /// - public void CutTrailZeros() + public void TrimTrailingZeros() { int count = 0; for (int ii = Degree - 1; ii >= 0; ii--) @@ -183,28 +185,28 @@ namespace MathNet.Numerics } var t = this.Clone() as Polynomial; - t.CutTrailZeros(); + t.TrimTrailingZeros(); var cNew = new double[t.Coeffs.Length - 1]; for (int i = 1; i < t.Coeffs.Length; i++) { cNew[i-1] = t.Coeffs[i] * i; } var p = new Polynomial(cNew, isFlip: IsFlipped); - p.CutTrailZeros(); + p.TrimTrailingZeros(); return p; } public Polynomial Integrate() { var t = this.Clone() as Polynomial; - t.CutTrailZeros(); + t.TrimTrailingZeros(); var cNew = new double[t.Coeffs.Length + 1]; for (int i = 1; i < cNew.Length; i++) { cNew[i] = t.Coeffs[i-1] / i; } var p = new Polynomial(cNew, isFlip: IsFlipped); - p.CutTrailZeros(); + p.TrimTrailingZeros(); return p; } @@ -222,13 +224,13 @@ namespace MathNet.Numerics public static Polynomial operator *( Polynomial a, Polynomial b) { // do not cut trailing zeros, since it may corrupt the outcom, if the array is of form 1 + x^-1 + x^-2 + x^-3 - //a.CutTrailZeros(); - //b.CutTrailZeros(); + //a.TrimTrailingZeros(); + //b.TrimTrailingZeros(); double[] ret = conv(a.Coeffs, b.Coeffs); Polynomial ret_p = new Polynomial(ret); - //ret_p.CutTrailZeros(); + //ret_p.TrimTrailingZeros(); return (ret_p); @@ -310,7 +312,7 @@ namespace MathNet.Numerics } /// - /// Calculates the complex roots of the Polynomial in the same way as matlab does + /// Calculates the complex roots of the Polynomial by eigenvalue decomposition /// /// a vector of complex numbers with the roots public Complex[] GetRoots() @@ -337,13 +339,14 @@ namespace MathNet.Numerics } /// - /// get the eigenvalue matrix A of this Polynomial such that eig(A) = roots of this Polynomial + /// get the eigenvalue matrix A of this Polynomial such that eig(A) = roots of this Polynomial. /// /// Eigenvalue matrix A + /// this matrix is similar to the companion matrix of this polynomial, in such a way, that it's transpose is the columnflip of the companion matrix public DenseMatrix GetEigValMatrix() { Polynomial pLoc = new Polynomial(this.Coeffs); - pLoc.CutTrailZeros(); + pLoc.TrimTrailingZeros(); int n = pLoc.Coeffs.Length - 1; if (n < 2) @@ -459,6 +462,112 @@ namespace MathNet.Numerics return (res_poly); } + /// + /// 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 + public static Tuple DivideLong(Polynomial a, Polynomial b) + { + if (a == null) + throw new ArgumentNullException("a"); + if (b == null) + throw new ArgumentNullException("b"); + + if (a.Degree <= 0) + throw new ArgumentOutOfRangeException("a Degree must be greater than zero"); + if (b.Degree <= 0) + throw new ArgumentOutOfRangeException("b Degree must be greater than zero"); + + if (b.Coeffs.Last() == 0) + throw new DivideByZeroException("b polynomial ends with zero"); + + var c1 = a.Coeffs.ToArray(); + var c2 = b.Coeffs.ToArray(); + + var n1 = c1.Length; + var n2 = c2.Length; + + double[] quo = null; + double[] rem = null; + + if (n2 == 1) // division by scalar + { + var fact = c2[0]; + quo = new double[n1]; + for (int i = 0; i < n1; i++) + quo[i] = c1[i] / fact; + rem = new double[] { 0 }; + } + else if(n1 < n2) // denominator degree higher than nominator degree + { + // quotient always be 0 and return c1 as remainder + quo = new double[] { 0 }; + rem = c1.ToArray(); + } + else + { + var dn = n1 - n2; + var scl = c2[n2 - 1]; + var c22 = new double[n2 - 1]; + for (int ii = 0; ii < c22.Length; ii++) + c22[ii] = c2[ii] / scl; + + int i = dn; + int j = n1 - 1; + while (i >= 0) + { + var vals = new double[j - i]; + for (int k = 0; k < vals.Length; k++) + vals[k] = c22[k] * c1[j]; + + for (int idx = i; idx < j; idx++) + c1[idx] -= vals[idx - i]; + + i--; + j--; + } + + rem = new double[j + 1]; + quo = new double[n1 - j + 1]; + + Array.Copy(c1, rem, j + 1); + + for (int idx2 = j+1; idx2 < n1; idx2++) + quo[idx2 - j + 1] = c1[idx2] / scl; + + } + + if (rem == null) + throw new NullReferenceException("resulting remainder was null"); + + if (quo == null) + throw new NullReferenceException("resulting quotient was null"); + + + // output mapping + var pQuo = new Polynomial(quo); + var pRem = new Polynomial(rem); + pQuo.TrimTrailingZeros(); + pQuo.TrimTrailingZeros(); + + return new Tuple(pQuo, pRem); + } + + + /// + /// 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 + public Tuple DivideLong(Polynomial b) + { + return DivideLong(this, b); + } + + #endregion #region Displaying @@ -492,7 +601,9 @@ namespace MathNet.Numerics for (int ii = 0; ii < Coeffs.Length; ii++) { - if (ii == 0) + if (ii == 0 && Coeffs.Length == 1) + strLoc += String.Format("{0}", this.Coeffs[ii], VarName, ii); + else if(ii == 0) strLoc += String.Format("{0} + ", this.Coeffs[ii], VarName, ii); else if (ii == Coeffs.Length - 1) strLoc += String.Format("{0}{1}{2}", this.Coeffs[ii], VarName, ii);