ufkc7e3d
/
mathnet-numerics
forked from tsai/mathnet-numerics
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

Polynomial: more refactoring

ridge-regression
Christoph Ruegg 8 years ago
parent
610921f6c2
commit
1ea8b9f32b
2 changed files with 320 additions and 280 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 22
      src/Numerics.Tests/PolynomialTests.cs
  2. 578
      src/Numerics/Polynomial.cs

22
src/Numerics.Tests/PolynomialTests.cs
View File

@ -228,25 +228,25 @@ namespace MathNet.Numerics.UnitTests
{ {
var p1 = new Polynomial(1.0d); var p1 = new Polynomial(1.0d);
var p2 = new Polynomial(new double[0]); var p2 = new Polynomial(new double[0]);
var tpl = Polynomial.DivideLong(p1, p2); var tpl = Polynomial.DivideRemainder(p1, p2);
}); });
Assert.Throws(typeof(ArgumentOutOfRangeException), () => Assert.Throws(typeof(ArgumentOutOfRangeException), () =>
{ {
var p1 = new Polynomial(1.0d); var p1 = new Polynomial(1.0d);
var p2 = new Polynomial(new double[0]); var p2 = new Polynomial(new double[0]);
var tpl = Polynomial.DivideLong(p2, p1); var tpl = Polynomial.DivideRemainder(p2, p1);
}); });
Assert.Throws(typeof(ArgumentOutOfRangeException), () => Assert.Throws(typeof(ArgumentOutOfRangeException), () =>
{ {
var p1 = new Polynomial(new double[0]); var p1 = new Polynomial(new double[0]);
var p2 = new Polynomial(new double[0]); var p2 = new Polynomial(new double[0]);
var tpl = Polynomial.DivideLong(p2, p1); var tpl = Polynomial.DivideRemainder(p2, p1);
}); });
Assert.Throws(typeof(DivideByZeroException), () => Assert.Throws(typeof(DivideByZeroException), () =>
{ {
var p1 = new Polynomial(1.0d); var p1 = new Polynomial(1.0d);
var p2 = new Polynomial(0.0d); var p2 = new Polynomial(0.0d);
var tpl = Polynomial.DivideLong(p1, p2); var tpl = Polynomial.DivideRemainder(p1, p2);
}); });
} }
@ -255,13 +255,13 @@ namespace MathNet.Numerics.UnitTests
{ {
var p11 = new Polynomial(2.0d); var p11 = new Polynomial(2.0d);
var p21 = new Polynomial(2.0d); var p21 = new Polynomial(2.0d);
var tpl1 = Polynomial.DivideLong(p11, p21); var tpl1 = Polynomial.DivideRemainder(p11, p21);
TestEqual(new double[] { 1.0 }, tpl1.Item1); TestEqual(new double[] { 1.0 }, tpl1.Item1);
TestEqual(new double[] { 0.0 }, tpl1.Item2); TestEqual(new double[] { 0.0 }, tpl1.Item2);
var p12 = new Polynomial(new double[] { 2.0d, 2.0d }); var p12 = new Polynomial(new double[] { 2.0d, 2.0d });
var p22 = new Polynomial(2.0d); var p22 = new Polynomial(2.0d);
var tpl2 = Polynomial.DivideLong(p12, p22); var tpl2 = Polynomial.DivideRemainder(p12, p22);
TestEqual(new double[] { 1.0, 1.0 }, tpl2.Item1); TestEqual(new double[] { 1.0, 1.0 }, tpl2.Item1);
TestEqual(new double[] { 0.0 }, tpl2.Item2); TestEqual(new double[] { 0.0 }, tpl2.Item2);
@ -280,7 +280,7 @@ namespace MathNet.Numerics.UnitTests
var pi = new Polynomial(ci); var pi = new Polynomial(ci);
var pj = new Polynomial(cj); var pj = new Polynomial(cj);
var tgt = Polynomial.Add(pi, pj); var tgt = Polynomial.Add(pi, pj);
var tpl3 = Polynomial.DivideLong(tgt, pi); var tpl3 = Polynomial.DivideRemainder(tgt, pi);
var pquo = tpl3.Item1; var pquo = tpl3.Item1;
var prem = tpl3.Item2; var prem = tpl3.Item2;
var pres = (pquo * pi) + prem; var pres = (pquo * pi) + prem;
@ -295,14 +295,14 @@ namespace MathNet.Numerics.UnitTests
public void GetRootsTest() public void GetRootsTest()
{ {
var tol = 1e-14; var tol = 1e-14;
// 0 = 1 -> no roots
var p1 = new Polynomial(1.0); var p1 = new Polynomial(1.0);
var r = p1.Roots(); var r = p1.Roots();
Assert.AreEqual(0, r.Length, "length mismatch");
Assert.AreEqual(1, r.Length, "length mismatch"); // 0 = 1 + 2*x -> single root at -1/2
Assert.AreEqual(1.0, r.FirstOrDefault().Real);
var p2 = new Polynomial(new double[] { 1, 2 }); var p2 = new Polynomial(new double[] { 1, 2 });
var r2 = p2.Roots(); var r2 = p2.Roots();
Assert.AreEqual(1, r2.Length, "length mismatch"); Assert.AreEqual(1, r2.Length, "length mismatch");
Assert.AreEqual(-0.5, r2.FirstOrDefault().Real, tol); Assert.AreEqual(-0.5, r2.FirstOrDefault().Real, tol);

578
src/Numerics/Polynomial.cs
View File

@ -25,78 +25,75 @@ namespace MathNet.Numerics
public string VarName = "x^"; public string VarName = "x^";
/// <summary> /// <summary>
/// Degree of the polynomial, i.e. the largest monomial exponent. For example, the degree of x^2+x^5 is 5. /// Degree of the polynomial, i.e. the largest monomial exponent. For example, the degree of y=x^2+x^5 is 5, for y=3 it is 0.
/// The null-polynomial returns degree -1 because the correct degree, negative infinity, cannot be represented by integers. /// The null-polynomial returns degree -1 because the correct degree, negative infinity, cannot be represented by integers.
/// </summary> /// </summary>
public int Degree public int Degree => EvaluateDegree(Coefficients);
{
get
{
if (Coefficients == null)
{
return -1;
}
for (int i = Coefficients.Length - 1; i >= 0; i--)
{
if (Coefficients[i] != 0.0)
{
return i;
}
}
return -1;
}
}
/// <summary> /// <summary>
/// constructor setting a Polynomial of size n containing only zeros /// Create a zero-polynomial with a coefficient array of the given length.
/// </summary> /// </summary>
/// <param name="n">size of Polynomial</param> /// <param name="n">Length of the coefficient array</param>
public Polynomial(int n) public Polynomial(int n)
{ {
if (n < 0) if (n < 0)
{ {
throw new ArgumentOutOfRangeException("n must be postive"); throw new ArgumentOutOfRangeException(nameof(n), "n must be non-negative");
} }
Coefficients = new double[n]; Coefficients = new double[n];
} }
/// <summary> /// <summary>
/// make Polynomial: e.G 3.0 = 3.0 + 0 x^1 + 0 x^2 /// Create a zero-polynomial
/// </summary> /// </summary>
/// <param name="coefficient">just the "x^0" part</param> public Polynomial()
public Polynomial(double coefficient)
{ {
Coefficients = new double[1]; Coefficients = new double[0];
Coefficients[0] = coefficient;
} }
/// <summary> /// <summary>
/// make Polynomial: e.G new double[] {5, 0, 2} = "5 + 0 x^1 + 2 x^2" /// Create a constant polynomial.
/// Example: 3.0 -> "p : x -> 3.0"
/// </summary> /// </summary>
/// <param name="coefficients">Polynomial coefficients as enumerable</param> /// <param name="coefficient">just the "x^0" part</param>
public Polynomial(IEnumerable<double> coefficients) public Polynomial(double coefficient)
{ {
if (coefficients == null) Coefficients = coefficient == 0.0 ? new double[0] : new[] { coefficient };
{
throw new ArgumentNullException(nameof(coefficients));
}
Coefficients = coefficients.ToArray();
} }
/// <summary> /// <summary>
/// make Polynomial: e.G new double[] {5, 0, 2} = "5 + 0 x^1 + 2 x^2" /// Create a polynomial with the provided coefficients (in ascending order, where the index matches the exponent).
/// Example: {5, 0, 2} -> "p : x -> 5 + 0 x^1 + 2 x^2".
/// </summary> /// </summary>
/// <param name="coefficients">Polynomial coefficients as array</param> /// <param name="coefficients">Polynomial coefficients as array</param>
public Polynomial(double[] coefficients) public Polynomial(double[] coefficients)
{ {
if (coefficients == null) int degree = EvaluateDegree(coefficients);
Coefficients = new double[degree + 1];
Array.Copy(coefficients, Coefficients, Coefficients.Length);
}
/// <summary>
/// Create a polynomial with the provided coefficients (in ascending order, where the index matches the exponent).
/// Example: {5, 0, 2} -> "p : x -> 5 + 0 x^1 + 2 x^2".
/// </summary>
/// <param name="coefficients">Polynomial coefficients as enumerable</param>
public Polynomial(IEnumerable<double> coefficients) : this(coefficients.ToArray())
{
}
static int EvaluateDegree(double[] coefficients)
{
for (int i = coefficients.Length - 1; i >= 0; i--)
{ {
throw new ArgumentNullException(nameof(coefficients)); if (coefficients[i] != 0.0)
{
return i;
}
} }
Coefficients = new double[coefficients.Length];
Array.Copy(coefficients, Coefficients, coefficients.Length); return -1;
} }
/// <summary> /// <summary>
@ -140,6 +137,23 @@ namespace MathNet.Numerics
return new Polynomial(coefficients); return new Polynomial(coefficients);
} }
/// <summary>
/// This method returns the coefficients of the Polynomial as an array the "IsFlipped" property,
/// which is set during construction is taken into account automatically.
/// </summary>
/// <returns>The coefficients of the polynomial as an array</returns>
public double[] ToArray()
{
return Coefficients.ToArray();
}
public object Clone()
{
// TODO: this assumes the constructor does a copy
return new Polynomial(Coefficients);
}
#region Evaluation
/// <summary> /// <summary>
/// Evaluate a polynomial at point x. /// Evaluate a polynomial at point x.
/// </summary> /// </summary>
@ -175,7 +189,9 @@ namespace MathNet.Numerics
{ {
return z.Select(Evaluate); return z.Select(Evaluate);
} }
#endregion
#region Calculus
public Polynomial Differentiate() public Polynomial Differentiate()
{ {
if (Coefficients.Length == 0) if (Coefficients.Length == 0)
@ -188,7 +204,7 @@ namespace MathNet.Numerics
var cNew = new double[t.Coefficients.Length - 1]; var cNew = new double[t.Coefficients.Length - 1];
for (int i = 1; i < t.Coefficients.Length; i++) for (int i = 1; i < t.Coefficients.Length; i++)
{ {
cNew[i-1] = t.Coefficients[i] * i; cNew[i - 1] = t.Coefficients[i] * i;
} }
var p = new Polynomial(cNew); var p = new Polynomial(cNew);
@ -210,186 +226,48 @@ namespace MathNet.Numerics
p.Trim(); p.Trim();
return p; return p;
} }
#endregion
/// <summary> #region Linear Algebra
/// Addition of two Polynomials (piecewise)
/// </summary>
/// <param name="a">Left polynomial</param>
/// <param name="b">Right polynomial</param>
/// <returns>Resulting Polynomial</returns>
public static Polynomial operator +(Polynomial a, Polynomial b)
{
return Add(a, b);
}
/// <summary>
/// adds a scalar to a polynomial.
/// </summary>
/// <param name="a">Polynomial</param>
/// <param name="k">Scalar value</param>
/// <returns>Resulting Polynomial</returns>
public static Polynomial operator +(Polynomial a, double k)
{
return Add(a, k);
}
/// <summary>
/// adds a scalar to a polynomial.
/// </summary>
/// <param name="k">Scalar value</param>
/// <param name="a">Polynomial</param>
/// <returns>Resulting Polynomial</returns>
public static Polynomial operator +(double k, Polynomial a)
{
return Add(a, k);
}
/// <summary>
/// Subtraction of two polynomial.
/// </summary>
/// <param name="a">Left polynomial</param>
/// <param name="b">Right polynomial</param>
/// <returns>Resulting Polynomial</returns>
public static Polynomial operator -(Polynomial a, Polynomial b)
{
return Subtract(a, b);
}
/// <summary>
/// Subtracts a scalar from a polynomial.
/// </summary>
/// <param name="a">Polynomial</param>
/// <param name="k">Scalar value</param>
/// <returns>Resulting Polynomial</returns>
public static Polynomial operator -(Polynomial a, double k)
{
return Subtract(a, k);
}
/// <summary>
/// Subtracts a polynomial from a scalar.
/// </summary>
/// <param name="k">Scalar value</param>
/// <param name="a">Polynomial</param>
/// <returns>Resulting Polynomial</returns>
public static Polynomial operator -(double k, Polynomial a)
{
return Subtract(k, a);
}
/// <summary>
/// Negates a polynomial.
/// </summary>
/// <param name="a">Polynomial</param>
/// <returns>Resulting Polynomial</returns>
public static Polynomial operator -(Polynomial a)
{
return Negate(a);
}
/// <summary>
/// multiplies a Polynomial by a Polynomial using convolution [ASINCO.libs.subfun.conv(a.Coeffs, b.Coeffs)]
/// </summary>
/// <param name="a">Left polynomial</param>
/// <param name="b">Right polynomial</param>
/// <returns>resulting Polynomial</returns>
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 = Convolution(aa.Coefficients, bb.Coefficients);
Polynomial result = new Polynomial(ret);
//ret_p.Trim();
return result;
}
/// <summary>
/// multiplies a Polynomial by a scalar
/// </summary>
/// <param name="a">Polynomial</param>
/// <param name="k">Scalar value</param>
/// <returns>Resulting Polynomial</returns>
public static Polynomial operator *(Polynomial a, double k)
{
var aa = a.Clone() as Polynomial;
for (int ii = 0; ii < aa.Coefficients.Length; ii++)
aa.Coefficients[ii] *= k;
return aa;
}
/// <summary>
/// divide Polynomial by scalar value
/// </summary>
/// <param name="a">Polynomial</param>
/// <param name="k">Scalar value</param>
/// <returns>Resulting Polynomial</returns>
public static Polynomial operator /(Polynomial a, double k)
{
var aa = a.Clone() as Polynomial;
for (int ii = 0; ii < aa.Coefficients.Length; ii++)
aa.Coefficients[ii] /= k;
return aa;
}
/// <summary> /// <summary>
/// Calculates the complex roots of the Polynomial by eigenvalue decomposition /// Calculates the complex roots of the Polynomial by eigenvalue decomposition
/// </summary> /// </summary>
/// <returns>a vector of complex numbers with the roots</returns> /// <returns>a vector of complex numbers with the roots</returns>
public Complex[] Roots() public Complex[] Roots()
{ {
DenseMatrix A = EigenvalueMatrix(); switch (Degree)
Complex[] roots;
if (A == null)
{ {
if (Coefficients.Length < 2) case -1: // Zero-polynomial
{ case 0: // Non-zero constant: y = a0
var val = Coefficients.Length == 1 ? Coefficients[0] : Double.NaN; return new Complex[0];
roots = new Complex[] { val }; case 1: // Linear: y = a0 + a1*x
} return new[] { new Complex(-Coefficients[0] / Coefficients[1], 0) };
else
roots = new[] { new Complex(-Coefficients[0] / Coefficients[1], 0) };
}
else
{
Evd<double> eigen = A.Evd(Symmetricity.Asymmetric);
roots = eigen.EigenValues.ToArray();
} }
return roots; DenseMatrix A = EigenvalueMatrix();
Evd<double> eigen = A.Evd(Symmetricity.Asymmetric);
return eigen.EigenValues.AsArray();
} }
/// <summary> /// <summary>
/// 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.
/// </summary> /// </summary>
/// <returns>Eigenvalue matrix A</returns> /// <returns>Eigenvalue matrix A</returns>
/// <note>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</note> /// <note>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</note>
public DenseMatrix EigenvalueMatrix() public DenseMatrix EigenvalueMatrix()
{ {
Polynomial pLoc = new Polynomial(Coefficients); int n = Degree;
pLoc.Trim();
int n = pLoc.Coefficients.Length - 1;
if (n < 2) if (n < 2)
{ {
return null; return null;
} }
double a0 = pLoc.Coefficients[n]; // Negate, and normalize (scale such that the polynomial becomes monic)
double aN = Coefficients[n];
double[] p = new double[n]; double[] p = new double[n];
for (int ii = n - 1; ii >= 0; ii--) for (int ii = n - 1; ii >= 0; ii--)
{ {
p[ii] = -pLoc.Coefficients[ii] / a0; p[ii] = -Coefficients[ii] / aN;
} }
DenseMatrix A0 = DenseMatrix.CreateDiagonal(n - 1, n - 1, 1.0); DenseMatrix A0 = DenseMatrix.CreateDiagonal(n - 1, n - 1, 1.0);
@ -400,7 +278,9 @@ namespace MathNet.Numerics
return A; return A;
} }
#endregion
#region Arithmetic Operations
/// <summary> /// <summary>
/// Addition of two Polynomials (point-wise). /// Addition of two Polynomials (point-wise).
/// </summary> /// </summary>
@ -544,61 +424,77 @@ namespace MathNet.Numerics
} }
/// <summary> /// <summary>
/// Point-wise division of two Polynomials /// Multiplies a polynomial by a polynomial (convolution)
/// </summary> /// </summary>
/// <param name="a">Left Polynomial</param> /// <param name="a">Left polynomial</param>
/// <param name="b">Right Polynomial</param> /// <param name="b">Right polynomial</param>
/// <returns>Resulting Polynomial</returns> /// <returns>Resulting Polynomial</returns>
public static Polynomial PointwiseDivide(Polynomial a, Polynomial b) public static Polynomial Multiply(Polynomial a, Polynomial b)
{ {
var ac = a.Coefficients; var aa = a.Clone() as Polynomial;
var bc = b.Coefficients; 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();
var degree = a.Degree; double[] a1 = aa.Coefficients;
var result = new double[degree + 1]; double[] b1 = bb.Coefficients;
double[] ret = new double[a1.Length + b1.Length];
var commonLength = Math.Min(Math.Min(ac.Length, bc.Length), result.Length); for (int i = 0; i < a1.Length; i++)
for (int i = 0; i < commonLength; i++)
{ {
result[i] = ac[i] / bc[i]; for (int j = 0; j < b1.Length; j++)
{
ret[i + j] += a1[i] * b1[j];
}
} }
for (int i = commonLength; i < result.Length; i++) Polynomial result = new Polynomial(ret);
{
result[i] = ac[i] / 0.0;
}
return new Polynomial(result); //ret_p.Trim();
return result;
} }
/// <summary> /// <summary>
/// Point-wise multiplication of two Polynomials /// Scales a polynomial by a scalar
/// </summary> /// </summary>
/// <param name="a">Left Polynomial</param> /// <param name="a">Polynomial</param>
/// <param name="b">Right Polynomial</param> /// <param name="k">Scalar value</param>
/// <returns>Resulting Polynomial</returns> /// <returns>Resulting Polynomial</returns>
public static Polynomial PointwiseMultiply(Polynomial a, Polynomial b) public static Polynomial Multiply(Polynomial a, double k)
{ {
var ac = a.Coefficients; var aa = a.Clone() as Polynomial;
var bc = b.Coefficients;
var degree = Math.Min(a.Degree, b.Degree); for (int ii = 0; ii < aa.Coefficients.Length; ii++)
var result = new double[degree + 1]; aa.Coefficients[ii] *= k;
for (int i = 0; i < result.Length; i++)
{
result[i] = ac[i] * bc[i];
}
return new Polynomial(result); return aa;
} }
/// <summary> /// <summary>
/// Division of two polynomials returning the quotient-with-remainder of the two polynomials given /// Scales a polynomial by division by a scalar
/// </summary>
/// <param name="a">Polynomial</param>
/// <param name="k">Scalar value</param>
/// <returns>Resulting Polynomial</returns>
public static Polynomial Divide(Polynomial a, double k)
{
var aa = a.Clone() as Polynomial;
for (int ii = 0; ii < aa.Coefficients.Length; ii++)
aa.Coefficients[ii] /= k;
return aa;
}
/// <summary>
/// Euclidean long division of two polynomials, returning the quotient q and remainder r of the two polynomials a and b such that a = q*b + r
/// </summary> /// </summary>
/// <param name="a">Left polynomial</param> /// <param name="a">Left polynomial</param>
/// <param name="b">Right polynomial</param> /// <param name="b">Right polynomial</param>
/// <returns>A tuple holding quotient in first and remainder in second</returns> /// <returns>A tuple holding quotient in first and remainder in second</returns>
public static Tuple<Polynomial, Polynomial> DivideLong(Polynomial a, Polynomial b) public static Tuple<Polynomial, Polynomial> DivideRemainder(Polynomial a, Polynomial b)
{ {
if (a == null) if (a == null)
throw new ArgumentNullException(nameof(a)); throw new ArgumentNullException(nameof(a));
@ -630,7 +526,7 @@ namespace MathNet.Numerics
quo[i] = c1[i] / fact; quo[i] = c1[i] / fact;
rem = new double[] { 0 }; rem = new double[] { 0 };
} }
else if(n1 < n2) // denominator degree higher than nominator degree else if (n1 < n2) // denominator degree higher than nominator degree
{ {
// quotient always be 0 and return c1 as remainder // quotient always be 0 and return c1 as remainder
quo = new double[] { 0 }; quo = new double[] { 0 };
@ -652,7 +548,7 @@ namespace MathNet.Numerics
{ {
var v = c1[j]; var v = c1[j];
for (int k = i; k < j; k++) for (int k = i; k < j; k++)
c1[k] -= c22[k-i] * v; c1[k] -= c22[k - i] * v;
i--; i--;
j--; j--;
} }
@ -689,24 +585,201 @@ namespace MathNet.Numerics
pQuo.Trim(); pQuo.Trim();
return new Tuple<Polynomial, Polynomial>(pQuo, pRem); return new Tuple<Polynomial, Polynomial>(pQuo, pRem);
} }
#endregion
#region Arithmetic Pointwise Operations
/// <summary>
/// Point-wise division of two Polynomials
/// </summary>
/// <param name="a">Left Polynomial</param>
/// <param name="b">Right Polynomial</param>
/// <returns>Resulting Polynomial</returns>
public static Polynomial PointwiseDivide(Polynomial a, Polynomial b)
{
var ac = a.Coefficients;
var bc = b.Coefficients;
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++)
{
result[i] = ac[i] / bc[i];
}
for (int i = commonLength; i < result.Length; i++)
{
result[i] = ac[i] / 0.0;
}
return new Polynomial(result);
}
/// <summary>
/// Point-wise multiplication of two Polynomials
/// </summary>
/// <param name="a">Left Polynomial</param>
/// <param name="b">Right Polynomial</param>
/// <returns>Resulting Polynomial</returns>
public static Polynomial PointwiseMultiply(Polynomial a, Polynomial b)
{
var ac = a.Coefficients;
var bc = b.Coefficients;
var degree = Math.Min(a.Degree, b.Degree);
var result = new double[degree + 1];
for (int i = 0; i < result.Length; i++)
{
result[i] = ac[i] * bc[i];
}
return new Polynomial(result);
}
#endregion
#region Arithmetic Instance Methods (forwarders)
/// <summary> /// <summary>
/// Division of two polynomials returning the quotient-with-remainder of the two polynomials given /// Division of two polynomials returning the quotient-with-remainder of the two polynomials given
/// </summary> /// </summary>
/// <param name="b">Right polynomial</param> /// <param name="b">Right polynomial</param>
/// <returns>A tuple holding quotient in first and remainder in second</returns> /// <returns>A tuple holding quotient in first and remainder in second</returns>
public Tuple<Polynomial, Polynomial> DivideLong(Polynomial b) public Tuple<Polynomial, Polynomial> DivideRemainder(Polynomial b)
{
return DivideRemainder(this, b);
}
#endregion
#region Arithmetic Operator Overloads (forwarders)
/// <summary>
/// Addition of two Polynomials (piecewise)
/// </summary>
/// <param name="a">Left polynomial</param>
/// <param name="b">Right polynomial</param>
/// <returns>Resulting Polynomial</returns>
public static Polynomial operator +(Polynomial a, Polynomial b)
{
return Add(a, b);
}
/// <summary>
/// adds a scalar to a polynomial.
/// </summary>
/// <param name="a">Polynomial</param>
/// <param name="k">Scalar value</param>
/// <returns>Resulting Polynomial</returns>
public static Polynomial operator +(Polynomial a, double k)
{
return Add(a, k);
}
/// <summary>
/// adds a scalar to a polynomial.
/// </summary>
/// <param name="k">Scalar value</param>
/// <param name="a">Polynomial</param>
/// <returns>Resulting Polynomial</returns>
public static Polynomial operator +(double k, Polynomial a)
{
return Add(a, k);
}
/// <summary>
/// Subtraction of two polynomial.
/// </summary>
/// <param name="a">Left polynomial</param>
/// <param name="b">Right polynomial</param>
/// <returns>Resulting Polynomial</returns>
public static Polynomial operator -(Polynomial a, Polynomial b)
{
return Subtract(a, b);
}
/// <summary>
/// Subtracts a scalar from a polynomial.
/// </summary>
/// <param name="a">Polynomial</param>
/// <param name="k">Scalar value</param>
/// <returns>Resulting Polynomial</returns>
public static Polynomial operator -(Polynomial a, double k)
{
return Subtract(a, k);
}
/// <summary>
/// Subtracts a polynomial from a scalar.
/// </summary>
/// <param name="k">Scalar value</param>
/// <param name="a">Polynomial</param>
/// <returns>Resulting Polynomial</returns>
public static Polynomial operator -(double k, Polynomial a)
{ {
return DivideLong(this, b); return Subtract(k, a);
} }
/// <summary>
/// Negates a polynomial.
/// </summary>
/// <param name="a">Polynomial</param>
/// <returns>Resulting Polynomial</returns>
public static Polynomial operator -(Polynomial a)
{
return Negate(a);
}
/// <summary>
/// Multiplies a polynomial by a polynomial (convolution).
/// </summary>
/// <param name="a">Left polynomial</param>
/// <param name="b">Right polynomial</param>
/// <returns>resulting Polynomial</returns>
public static Polynomial operator *(Polynomial a, Polynomial b)
{
return Multiply(a, b);
}
/// <summary>
/// Multiplies a polynomial by a scalar.
/// </summary>
/// <param name="a">Polynomial</param>
/// <param name="k">Scalar value</param>
/// <returns>Resulting Polynomial</returns>
public static Polynomial operator *(Polynomial a, double k)
{
return Multiply(a, k);
}
/// <summary>
/// Multiplies a polynomial by a scalar.
/// </summary>
/// <param name="k">Scalar value</param>
/// <param name="a">Polynomial</param>
/// <returns>Resulting Polynomial</returns>
public static Polynomial operator *(double k, Polynomial a)
{
return Multiply(a, k);
}
/// <summary>
/// Divides a polynomial by scalar value.
/// </summary>
/// <param name="a">Polynomial</param>
/// <param name="k">Scalar value</param>
/// <returns>Resulting Polynomial</returns>
public static Polynomial operator /(Polynomial a, double k)
{
return Divide(a, k);
}
#endregion
#region ToString
/// <summary> /// <summary>
/// "0.00 x^3 + 0.00 x^2 + 0.00 x^1 + 0.00" like display of this Polynomial /// "0.00 x^3 + 0.00 x^2 + 0.00 x^1 + 0.00" like display of this Polynomial
/// </summary> /// </summary>
/// <returns>string in displayed format</returns> /// <returns>string in displayed format</returns>
public override string ToString() public override string ToString()
{ {
return ToString(highestFirst:false); return ToString(highestFirst: false);
} }
/// <summary> /// <summary>
@ -732,7 +805,7 @@ namespace MathNet.Numerics
if (ii == 0 && Coefficients.Length == 1) if (ii == 0 && Coefficients.Length == 1)
result += String.Format("{0}", Coefficients[ii], VarName, ii); result += String.Format("{0}", Coefficients[ii], VarName, ii);
else if(ii == 0) else if (ii == 0)
result += String.Format("{0} + ", Coefficients[ii], VarName, ii); result += String.Format("{0} + ", Coefficients[ii], VarName, ii);
else if (ii == Coefficients.Length - 1) else if (ii == Coefficients.Length - 1)
result += String.Format("{0}{1}{2}", Coefficients[ii], VarName, ii); result += String.Format("{0}{1}{2}", Coefficients[ii], VarName, ii);
@ -753,39 +826,6 @@ namespace MathNet.Numerics
return result; return result;
} }
#endregion
/// <summary>
/// This method returns the coefficients of the Polynomial as an array the "IsFlipped" property,
/// which is set during construction is taken into account automatically.
/// </summary>
/// <returns>the coefficients of the Polynomial as an array</returns>
public double[] ToArray()
{
return Coefficients.ToArray();
}
/// <summary>
/// Full convolution of two arrays
/// </summary>
/// <returns>convolution of a and b as vector</returns>
static double[] Convolution(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()
{
// TODO: this assumes the constructor does a copy
return new Polynomial(Coefficients);
}
} }
} }

Write Preview
Loading…
Cancel
Save