diff --git a/src/Numerics.Tests/PolynomialTests.cs b/src/Numerics.Tests/PolynomialTests.cs new file mode 100644 index 00000000..9b96f3be --- /dev/null +++ b/src/Numerics.Tests/PolynomialTests.cs @@ -0,0 +1,394 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// +// Copyright (c) 2009-2018 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using System; +using System.Collections.Generic; +using System.Diagnostics; +using System.Linq; +using System.Numerics; +using MathNet.Numerics; +using MathNet.Numerics.LinearRegression; +using MathNet.Numerics.Statistics; +using NUnit.Framework; + +namespace MathNet.Numerics.UnitTests +{ + /// + /// some of these tests were inspired by numpys tests in python for the Polynomial functions. + /// Thanks to the numpy contributers! + /// + [TestFixture, Category("Calculus")] + public class PolynomialTests + { + + [TestCase(new double[] { 5, 4, 3, 0, 2 }, "5 + 4x^1 + 3x^2 + 0x^3 + 2x^4")] + [TestCase(new double[0], "")] + [TestCase(new double[] { 0, 4, 3, 0, 0 }, "0 + 4x^1 + 3x^2 + 0x^3 + 0x^4")] + public void ToStringTest(double[] x, string expected) + { + var p = new Polynomial(x); + Assert.AreEqual(expected, p.ToString()); + } + + [TestCase(new double[] { 5, 4, 3, 0, 2 }, new double[] { 4*1, 3*2, 0*3, 2*4 })] + [TestCase(new double[0], null)] + [TestCase(new double[] { 0, 4, 3, 0, 0 }, new double[] { 4*1, 3*2 })] + public void DifferentiateTest(double[] x, double[] expected) + { + var p = new Polynomial(x); + var p_res = p.Differentiate(); + + if (expected == null) + { + Assert.IsNull(p_res); + return; + } + else + { + Assert.AreEqual(expected.Length, p_res.Degree, "length mismatch"); + for (int k = 0; k < p_res.Degree; k++) + { + Assert.AreEqual(expected[k], p_res.Coeffs[k], "idx: " + k + " mismatch"); + } + } + + } + + [TestCase(new double[] { 5, 4, 3, 0, 2 }, new double[] { 0, 5.0/1.0, 4.0/2.0, 3.0/3.0, 0.0/4.0, 2.0/5.0 })] + [TestCase(new double[0], new double[1] { 0 })] + [TestCase(new double[] { 0, 1, 6, 8 }, new double[] {0, 0.0/1.0, 1.0/2.0, 6.0/3.0, 8.0/4.0})] + public void IntegrateTest(double[] x, double[] expected) + { + var p = new Polynomial(x); + var p_res = p.Integrate(); + + if (expected == null) + { + Assert.IsNull(p_res); + return; + } + else + { + Assert.AreEqual(expected.Length, p_res.Degree, "length mismatch"); + for (int k = 0; k < p_res.Degree; k++) + { + Assert.AreEqual(expected[k], p_res.Coeffs[k], "idx: " + k + " mismatch"); + } + } + } + + [Test] + public void AddTest() + { + 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 n = Math.Max(i, j) + 1; + var tgt = new double[n]; + tgt[i] += 1; + tgt[j] += 1; + + var c1 = new double[i + 1]; + var c2 = new double[j + 1]; + + c1[i] = 1.0; + c2[j] = 1.0; + + var p1 = new Polynomial(c1); + var p2 = new Polynomial(c2); + + var p_res = Polynomial.Add(p1, p2); + var p_tar = new Polynomial(tgt); + + p_res.Trim(); + p_tar.Trim(); + + 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); + } + } + } + } + + [Test] + public void SubstractTest() + { + 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 n = Math.Max(i, j) + 1; + var tgt = new double[n]; + tgt[i] += 1; + tgt[j] -= 1; + + var c1 = new double[i + 1]; + var c2 = new double[j + 1]; + + c1[i] = 1.0; + c2[j] = 1.0; + + var p1 = new Polynomial(c1); + var p2 = new Polynomial(c2); + + var p_res = Polynomial.Substract(p1, p2); + var p_tar = new Polynomial(tgt); + + p_res.Trim(); + p_tar.Trim(); + + 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); + } + } + } + } + + [Test] + public void MultiplyTest() + { + 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 n = i + j + 1; + var tgt = new double[n]; + tgt[i + j] += 1; + + var c1 = new double[i + 1]; + var c2 = new double[j + 1]; + + c1[i] = 1.0; + c2[j] = 1.0; + + var p1 = new Polynomial(c1); + var p2 = new Polynomial(c2); + + var p_res = p1 * p2; + var p_tar = new Polynomial(tgt); + + p_res.Trim(); + p_tar.Trim(); + + 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); + } + } + } + } + + + [TestCase(new double[] { 5, 4, 0 }, "5 + 4x^1")] + [TestCase(new double[] { 0, 0, 0 }, "0")] + [TestCase(new double[] { 5, 4, 3, 0, 2 }, "5 + 4x^1 + 3x^2 + 0x^3 + 2x^4")] + [TestCase(new double[] { 0, 0, 8, 0, 0 }, "0 + 0x^1 + 8x^2")] + [TestCase(new double[] { 0, 4, 3, 0, 0 }, "0 + 4x^1 + 3x^2")] + public void TrimTest(double[] x, string expected) + { + var p = new Polynomial(x); + p.Trim(); + Assert.AreEqual(expected, p.ToString()); + + } + + 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 DivideLongTestWrongInputs() + { + 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); + }); + } + + [Test] + public void DivideLongTest() + { + + 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[i + 2]; + var cj = new double[j + 2]; + 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; + pres.Trim(); + + testEqual(pres, tgt, msg); + } + } + } + + [Test] + public void GetRootsTest() + { + var tol = 1e-14; + var p1 = new Polynomial(1.0); + var r = p1.GetRoots(); + + Assert.AreEqual(1, r.Length, "length mismatch"); + Assert.AreEqual(1.0, r.FirstOrDefault().Real); + + var p2 = new Polynomial(new double[] { 1, 2 }); + + var r2 = p2.GetRoots(); + Assert.AreEqual(1, r2.Length, "length mismatch"); + Assert.AreEqual(-0.5, r2.FirstOrDefault().Real, tol); + + // T.G: the following expected values were generated using + // numpys np.roots(x) method + // which is equivalent to np.polynomial.polynomial.polyroots + + 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); + + 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); + + 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); + } + + private void testEqual(double[] x, List eIn) + { + var tol = 1e-10; + var r0 = new Polynomial(x).GetRoots().ToList(); + + var e = eIn.OrderBy(v => v.Real).ToArray(); + var r = r0.OrderBy(v => v.Real).ToArray(); + + Assert.IsNotNull(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); + } + } + + 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/FindRoots.cs b/src/Numerics/FindRoots.cs index 3e8af9a7..d8dc8963 100644 --- a/src/Numerics/FindRoots.cs +++ b/src/Numerics/FindRoots.cs @@ -110,6 +110,15 @@ namespace MathNet.Numerics { return RootFinding.Cubic.Roots(d, c, b, a); } + /// + /// Find all roots of a polynomial by calculating the characteristic polynomial of the companion matrix + /// + /// the values for the polynomial in ascending order e.G new double[] {5, 0, 2} = "5 + 0 x^1 + 2 x^2" + /// the roots of the polynomial + public static Complex[] Polynomial(double[] poly) + { + return new Polynomial(poly).GetRoots(); + } /// /// Find all roots of the Chebychev polynomial of the first kind. diff --git a/src/Numerics/Polynomial.cs b/src/Numerics/Polynomial.cs new file mode 100644 index 00000000..64a57887 --- /dev/null +++ b/src/Numerics/Polynomial.cs @@ -0,0 +1,696 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using System.Numerics; +using MathNet.Numerics; +using MathNet.Numerics.LinearAlgebra; +using MathNet.Numerics.LinearAlgebra.Double; +using MathNet.Numerics.Statistics; +using MathNet.Numerics.IntegralTransforms; +using MathNet.Numerics.LinearRegression; +using MathNet.Numerics.LinearAlgebra.Factorization; + +namespace MathNet.Numerics +{ + /// + /// 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; } + + /// + /// 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 Degree + { + get + { + return (Coeffs == null ? 0 : Coeffs.Length); + } + } + + /// + /// constructor setting a Polynomial of size n containing only zeros + /// + /// size of Polynomial + public Polynomial(int n) + { + if (n < 0) + { + throw new ArgumentOutOfRangeException("n must be postive"); + } + Coeffs = new double[n]; + } + + + /// + /// make Polynomial: e.G 3.0 = 3.0 + 0 x^1 + 0 x^2 + /// + /// just the "x^0" part + public Polynomial(double coeff) + { + this.Coeffs = new double[1]; + Coeffs[0] = coeff; + } + + /// + /// make Polynomial: e.G new double[] {5, 0, 2} = "5 + 0 x^1 + 2 x^2" + /// + /// Polynomial coefficiens as enumerable + public Polynomial(IEnumerable coeffs) + { + if (coeffs == null) + { + throw new ArgumentNullException("coeffs"); + } + this.Coeffs = coeffs.ToArray(); + } + + /// + /// make Polynomial: e.G new double[] {5, 0, 2} = "5 + 0 x^1 + 2 x^2" + /// + /// Polynomial coefficiens as array + public Polynomial(double[] coeffs) + { + if (coeffs == null) + { + throw new ArgumentNullException("coeffs"); + } + this.Coeffs = new double[coeffs.Length]; + Array.Copy(coeffs, this.Coeffs, coeffs.Length); + } + + /// + /// remove all trailing zeros, e.G before: "3 + 2 x^1 + 0 x^2" after: "3 + 2 x^1" + /// + public void Trim() + { + if (Degree == 1) + return; + + int i = Degree - 1; + while (i >= 0 && Coeffs[i] == 0.0) + i--; + + if (i < 0) + Coeffs = new double[1] { 0.0 }; + else if (i == 0) + Coeffs = new double[1] { Coeffs[0] }; + else + { + var hold = new double[i+1]; + Array.Copy(Coeffs, hold, i+1); + Coeffs = hold; + } + } + + + #region Data Interaction + + /// + /// Least-Squares fitting the points (x,y) to a k-order polynomial y : x -> p0 + p1*x + p2*x^2 + ... + pk*x^k + /// + public static Polynomial Fit(double[] x, double[] y, int order, DirectRegressionMethod method = DirectRegressionMethod.QR) + { + var pArr = MathNet.Numerics.Fit.Polynomial(x, y, order, method); + return new Polynomial(pArr); + } + + /// + /// Evaluate a polynomial at point x. + /// + /// The location where to evaluate the polynomial at. + public double Evaluate(double z) + { + return MathNet.Numerics.Evaluate.Polynomial(z, Coeffs); + } + + /// + /// Evaluate a polynomial at points z. + /// + /// The locations where to evaluate the polynomial at. + public IEnumerable Evaluate(IEnumerable z) + { + var Lst = new List(); + foreach (var item in z) + { + Lst.Add(Evaluate(item)); + } + return Lst; + } + + #endregion + + #region diff/int + public Polynomial Differentiate() + { + + if (Coeffs.Length == 0) + { + return null; + } + + var t = this.Clone() as Polynomial; + t.Trim(); + 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); + p.Trim(); + return p; + } + + public Polynomial Integrate() + { + var t = this.Clone() as Polynomial; + t.Trim(); + 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); + p.Trim(); + return p; + } + + #endregion + + #region Operators + + + /// + /// multiplies a Polynomial by a Polynomial using convolution [ASINCO.libs.subfun.conv(a.Coeffs, b.Coeffs)] + /// + /// left Polynomial + /// right Polynomial + /// resulting Polynomial + public static Polynomial operator *( Polynomial a, Polynomial b) + { + var aa = a.Clone() as Polynomial; + var bb = b.Clone() as Polynomial; + // 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.Trim(); + //b.Trim(); + + double[] ret = conv(aa.Coeffs, bb.Coeffs); + Polynomial ret_p = new Polynomial(ret); + + //ret_p.Trim(); + + return (ret_p); + + } + + /// + /// multiplies a Polynomial by a scalar + /// + /// left Polynomial + /// scalar value + /// resulting Polynomial + public static Polynomial operator *( Polynomial a, double k) + { + var aa = a.Clone() as Polynomial; + + + for (int ii = 0; ii < aa.Coeffs.Length; ii++) + aa.Coeffs[ii] *= k; + + return aa; + } + + /// + /// adds a scalar to a Polynomial (to the x^0 element) + /// + /// left Polynomial + /// scalar value + /// resulting Polynomial + public static Polynomial operator +( Polynomial a, double k) + { + var aa = a.Clone() as Polynomial; + + aa.Coeffs[0] += k; + return aa; + } + + /// + /// substracs a scalar from a Polynomial (from the x^0 element) + /// + /// left Polynomial + /// scalar value + /// resulting Polynomial + public static Polynomial operator -( Polynomial a, double k) + { + var aa = a.Clone() as Polynomial; + + a.Coeffs[0] -= k; + return aa; + } + + /// + /// divide Polynomial by scalar value + /// + /// left Polynomial + /// scalar value + /// resulting Polynomial + public static Polynomial operator /( Polynomial a, double k) + { + var aa = a.Clone() as Polynomial; + + for (int ii = 0; ii < aa.Coeffs.Length; ii++) + aa.Coeffs[ii] /= k; + + return aa; + } + + /// + /// Addition of two Polynomials (piecewise) + /// + /// left Polynomial + /// right Polynomial + /// resulting Polynomial + public static Polynomial operator +( Polynomial a, Polynomial b) + { + return Add(a, b); + } + + /// + /// substraction of two Polynomials (piecewise) + /// + /// left Polynomial + /// right Polynomial + /// resulting Polynomial + public static Polynomial operator -( Polynomial a, Polynomial b) + { + return Substract(a, b); + } + + /// + /// Calculates the complex roots of the Polynomial by eigenvalue decomposition + /// + /// a vector of complex numbers with the roots + public Complex[] GetRoots() + { + DenseMatrix A = this.GetEigValMatrix(); + Complex[] c_vec; + + if (A == null) + { + if (Coeffs.Length < 2) + { + var val = Coeffs.Length == 1 ? Coeffs[0] : Double.NaN; + c_vec = new Complex[1] { val }; + } + else + c_vec = new Complex[1] { new Complex(-Coeffs[0] / Coeffs[1], 0) }; + } + else + { + Evd eigen = A.Evd(Symmetricity.Asymmetric); + c_vec = eigen.EigenValues.ToArray(); + } + return c_vec; + } + + /// + /// 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.Trim(); + + int n = pLoc.Coeffs.Length - 1; + if (n < 2) + return null; + + double[] p = new double[n]; + + double a0 = pLoc.Coeffs[n]; + + for (int ii = n - 1; ii >= 0; ii--) + p[ii] = -pLoc.Coeffs[ii] / a0; + + 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()); + return A; + } + + /// + /// pointwise division of two Polynomials + /// + /// left Polynomial + /// right Polynomial + /// resulting Polynomial + public static Polynomial DividePointwise( Polynomial a, Polynomial b) + { + var aa = a.Clone() as Polynomial; + var bb = b.Clone() as Polynomial; + + if (aa.Coeffs.Length != bb.Coeffs.Length) + mkSameLength(ref aa, ref bb); + + int n = aa.Coeffs.Length; + double[] res = new double[aa.Coeffs.Length]; + + + for (int ii = 0; ii < n; ii++) + { + res[ii] = aa.Coeffs[ii] / bb.Coeffs[ii]; + } + Polynomial res_poly = new Polynomial(res); + return (res_poly); + } + + /// + /// pointwise multiplication of two Polynomials + /// + /// left Polynomial + /// right Polynomial + /// resulting Polynomial + public static Polynomial MultiplyPointwise( Polynomial a, Polynomial b) + { + var aa = a.Clone() as Polynomial; + var bb = b.Clone() as Polynomial; + + if (aa.Coeffs.Length != bb.Coeffs.Length) + mkSameLength(ref aa, ref bb); + + int n = aa.Coeffs.Length; + double[] res = new double[aa.Coeffs.Length]; + + + for (int ii = 0; ii < n; ii++) + { + res[ii] = aa.Coeffs[ii] * bb.Coeffs[ii]; + } + Polynomial res_poly = new Polynomial(res); + return (res_poly); + } + + /// + /// Addition of two Polynomials (piecewise) + /// + /// 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; + + if (aa.Degree != bb.Degree) + mkSameLength(ref aa, ref bb); + + int n = aa.Degree; + double[] res = new double[n]; + + + for (int ii = 0; ii < n; ii++) + { + res[ii] = aa.Coeffs[ii] + bb.Coeffs[ii]; + } + Polynomial res_poly = new Polynomial(res); + return (res_poly); + } + + /// + /// substraction of two Polynomials (piecewise) + /// + /// left Polynomial + /// right Polynomial + /// resulting Polynomial + public static Polynomial Substract( Polynomial a, Polynomial b) + { + var aa = a.Clone() as Polynomial; + var bb = b.Clone() as Polynomial; + + if (aa.Degree != bb.Degree) + mkSameLength(ref aa, ref bb); + + int n = aa.Degree; + double[] res = new double[n]; + + + for (int ii = 0; ii < n; ii++) + { + res[ii] = aa.Coeffs[ii] - bb.Coeffs[ii]; + } + Polynomial res_poly = new Polynomial(res); + 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[b.Degree-1] == 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 v = c1[j]; + var vals = new double[j - i]; + for (int k = i; k < j; k++) + c1[k] -= c22[k-i] * v; + i--; + j--; + } + + var j1 = j + 1; + var l1 = n1 - j1; + + rem = new double[j1]; + quo = new double[l1]; + + for (int k = 0; k < l1; k++) + quo[k] = c1[k + j1] / scl; + + for (int k = 0; k < j1; k++) + rem[k] = c1[k]; + + } + + 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); + + pRem.Trim(); + pQuo.Trim(); + return new Tuple(pQuo, pRem); + } + + + /// + /// 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 + public Tuple DivideLong(Polynomial b) + { + return DivideLong(this, b); + } + + + #endregion + + #region Displaying + /// + /// "0.00 x^3 + 0.00 x^2 + 0.00 x^1 + 0.00" like display of this Polynomial + /// + /// string in displayed format + public override string ToString() + { + return ToString(highestFirst:false); + } + + /// + /// "0.00 x^3 + 0.00 x^2 + 0.00 x^1 + 0.00" like display of this Polynomial + /// + /// string in displayed format + public string ToString(bool highestFirst) + { + string strLoc = ""; + if (this.Coeffs == null) + { + return "null"; + } + if (this.Coeffs.Length == 0) + { + return ""; + } + + if (!highestFirst) + { + for (int ii = 0; ii < Coeffs.Length; ii++) + { + + 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); + else + strLoc += String.Format("{0}{1}{2} + ", this.Coeffs[ii], VarName, ii); + } + } + else + { + for (int ii = Coeffs.Length - 1; ii >= 0; ii--) + { + if (ii == 0) + strLoc += this.Coeffs[ii].ToString(); + else + strLoc += String.Format("{0}{1}{2} + ", this.Coeffs[ii], VarName, ii); + } + } + + return strLoc; + } + + + #endregion + + #region Interfacing + + /// + /// This method returns the coefficcients of the Polynomial as an array the "IsFlipped" property, + /// which is set during construction is taken into account automatically. + /// + /// the coefficcients of the Polynomial as an array + public double[] ToArray() + { + return (Coeffs.ToArray()); + } + + #endregion + + #region Helpers + + private static void mkSameLength(ref Polynomial a, ref Polynomial b) + { + double[] aHold = new double[a.Coeffs.Length]; + double[] bHold = new double[b.Coeffs.Length]; + Array.Copy(a.Coeffs, aHold, a.Coeffs.Length); + Array.Copy(b.Coeffs, bHold, b.Coeffs.Length); + + if (a.Coeffs.Length < b.Coeffs.Length) + { + a.Coeffs = new double[b.Coeffs.Length]; + b.Coeffs = new double[b.Coeffs.Length]; + Array.Copy(aHold, a.Coeffs, aHold.Length); + Array.Copy(bHold, b.Coeffs, bHold.Length); + } + else + { + a.Coeffs = new double[a.Coeffs.Length]; + b.Coeffs = new double[a.Coeffs.Length]; + Array.Copy(aHold, a.Coeffs, aHold.Length); + Array.Copy(bHold, b.Coeffs, bHold.Length); + } + + } + + /// + /// (full) convolution of two arrays + /// + /// left vector + /// right vector + /// convolution of a and b as vector + private static double[] conv(double[] a, double[] b) + { + double[] ret = new double[a.Length + b.Length]; + + for (int i = 0; i < a.Length; i++) + { + for (int j = 0; j < b.Length; j++) + { + ret[i + j] += a[i] * b[j]; + } + } + return ret; + } + + public object Clone() + { + return new Polynomial(this.Coeffs); + } + #endregion + + } + +}