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

Polynomial: reimplement arithmetics in immutable way

ridge-regression
Christoph Ruegg 8 years ago
parent
674a221956
commit
9a641040f9
2 changed files with 128 additions and 135 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. 57
      src/Numerics.Tests/PolynomialTests.cs
  2. 206
      src/Numerics/Polynomial.cs

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

@ -66,8 +66,8 @@ namespace MathNet.Numerics.UnitTests
} }
else else
{ {
Assert.AreEqual(expected.Length, p_res.CoefficientCount, "length mismatch"); Assert.AreEqual(expected.Length, p_res.Coefficients.Length, "length mismatch");
for (int k = 0; k < p_res.CoefficientCount; k++) for (int k = 0; k < p_res.Coefficients.Length; k++)
{ {
Assert.AreEqual(expected[k], p_res.Coefficients[k], "idx: " + k + " mismatch"); Assert.AreEqual(expected[k], p_res.Coefficients[k], "idx: " + k + " mismatch");
} }
@ -89,8 +89,8 @@ namespace MathNet.Numerics.UnitTests
} }
else else
{ {
Assert.AreEqual(expected.Length, p_res.CoefficientCount, "length mismatch"); Assert.AreEqual(expected.Length, p_res.Coefficients.Length, "length mismatch");
for (int k = 0; k < p_res.CoefficientCount; k++) for (int k = 0; k < p_res.Coefficients.Length; k++)
{ {
Assert.AreEqual(expected[k], p_res.Coefficients[k], "idx: " + k + " mismatch"); Assert.AreEqual(expected[k], p_res.Coefficients[k], "idx: " + k + " mismatch");
} }
@ -125,8 +125,8 @@ namespace MathNet.Numerics.UnitTests
p_res.Trim(); p_res.Trim();
p_tar.Trim(); p_tar.Trim();
Assert.AreEqual(p_tar.CoefficientCount, p_res.CoefficientCount, "length mismatch"); Assert.AreEqual(p_tar.Coefficients.Length, p_res.Coefficients.Length, "length mismatch");
for (int k = 0; k < p_res.CoefficientCount; k++) for (int k = 0; k < p_res.Coefficients.Length; k++)
{ {
Assert.AreEqual(p_tar.Coefficients[k], p_res.Coefficients[k], msg); Assert.AreEqual(p_tar.Coefficients[k], p_res.Coefficients[k], msg);
} }
@ -135,7 +135,7 @@ namespace MathNet.Numerics.UnitTests
} }
[Test] [Test]
public void SubstractTest() public void SubtractTest()
{ {
for (int i = 0; i < 5; i++) for (int i = 0; i < 5; i++)
{ {
@ -156,14 +156,14 @@ namespace MathNet.Numerics.UnitTests
var p1 = new Polynomial(c1); var p1 = new Polynomial(c1);
var p2 = new Polynomial(c2); 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); var p_tar = new Polynomial(tgt);
p_res.Trim(); p_res.Trim();
p_tar.Trim(); p_tar.Trim();
Assert.AreEqual(p_tar.CoefficientCount, p_res.CoefficientCount, "length mismatch"); Assert.AreEqual(p_tar.Coefficients.Length, p_res.Coefficients.Length, "length mismatch");
for (int k = 0; k < p_res.CoefficientCount; k++) for (int k = 0; k < p_res.Coefficients.Length; k++)
{ {
Assert.AreEqual(p_tar.Coefficients[k], p_res.Coefficients[k], msg); Assert.AreEqual(p_tar.Coefficients[k], p_res.Coefficients[k], msg);
} }
@ -198,8 +198,8 @@ namespace MathNet.Numerics.UnitTests
p_res.Trim(); p_res.Trim();
p_tar.Trim(); p_tar.Trim();
Assert.AreEqual(p_tar.CoefficientCount, p_res.CoefficientCount, "length mismatch"); Assert.AreEqual(p_tar.Coefficients.Length, p_res.Coefficients.Length, "length mismatch");
for (int k = 0; k < p_res.CoefficientCount; k++) for (int k = 0; k < p_res.Coefficients.Length; k++)
{ {
Assert.AreEqual(p_tar.Coefficients[k], p_res.Coefficients[k], msg); 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 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.DivideLong(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.DivideLong(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);
for (int i = 0; i < 5; i++) for (int i = 0; i < 5; i++)
{ {
@ -286,7 +286,7 @@ namespace MathNet.Numerics.UnitTests
var pres = (pquo * pi) + prem; var pres = (pquo * pi) + prem;
pres.Trim(); 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 x_2 = new double[] { -1.0, 1.0 };
var expected_2 = new List<Complex>(); var expected_2 = new List<Complex>();
expected_2.Add(new Complex(1.0, 0.0)); 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 x_3 = new double[] { -1.0, 0.0, 1.0 };
var expected_3 = new List<Complex>(); var expected_3 = new List<Complex>();
expected_3.Add(new Complex(1.0, 0.0)); expected_3.Add(new Complex(1.0, 0.0));
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 x_4 = new double[] { -1.0, -0.33333333333333337, 0.33333333333333326, 1.0 };
var expected_4 = new List<Complex>(); var expected_4 = new List<Complex>();
expected_4.Add(new Complex(0.9999999999999996, 0.0)); 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));
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<Complex> eIn) static void TestEqual(double[] x, List<Complex> eIn)
{ {
var tol = 1e-10; var tol = 1e-10;
var r0 = new Polynomial(x).GetRoots().ToList(); 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"); Assert.AreEqual(p_tar.Length, p_res.Length, "length mismatch");
for (int k = 0; k < p_res.Length; k++) 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"); Assert.AreEqual(p_tar.Length, p_res.Coefficients.Length, "length mismatch");
for (int k = 0; k < p_res.CoefficientCount; k++) for (int k = 0; k < p_res.Coefficients.Length; k++)
{ {
Assert.AreEqual(p_tar[k], p_res.Coefficients[k], msg); 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"); Assert.AreEqual(p_tar.Coefficients.Length, p_res.Coefficients.Length, "length mismatch");
for (int k = 0; k < p_res.CoefficientCount; k++) for (int k = 0; k < p_res.Coefficients.Length; k++)
{ {
Assert.AreEqual(p_tar.Coefficients[k], p_res.Coefficients[k], msg); Assert.AreEqual(p_tar.Coefficients[k], p_res.Coefficients[k], msg);
} }

206
src/Numerics/Polynomial.cs
View File

@ -17,24 +17,13 @@ namespace MathNet.Numerics
/// <summary> /// <summary>
/// The coefficients of the polynomial in a /// The coefficients of the polynomial in a
/// </summary> /// </summary>
public double[] Coefficients { get; set; } public double[] Coefficients { get; private set; }
/// <summary> /// <summary>
/// Only needed for the ToString method /// Only needed for the ToString method
/// </summary> /// </summary>
public string VarName = "x^"; public string VarName = "x^";
/// <summary>
/// Length of Polynomial (max element + 1) e.G x^5 highest element, will give Length = 6
/// </summary>
public int CoefficientCount
{
get
{
return (Coefficients == null ? 0 : Coefficients.Length);
}
}
/// <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 x^2+x^5 is 5.
/// 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.
@ -115,17 +104,25 @@ namespace MathNet.Numerics
/// </summary> /// </summary>
public void Trim() public void Trim()
{ {
if (CoefficientCount == 1) if (Coefficients.Length == 1)
{
return; return;
}
int i = CoefficientCount - 1; int i = Coefficients.Length - 1;
while (i >= 0 && Coefficients[i] == 0.0) while (i >= 0 && Coefficients[i] == 0.0)
{
i--; i--;
}
if (i < 0) if (i < 0)
{
Coefficients = new[] { 0.0 }; Coefficients = new[] { 0.0 };
}
else if (i == 0) else if (i == 0)
{
Coefficients = new[] { Coefficients[0] }; Coefficients = new[] { Coefficients[0] };
}
else else
{ {
var hold = new double[i+1]; var hold = new double[i+1];
@ -304,7 +301,7 @@ namespace MathNet.Numerics
/// <returns>resulting Polynomial</returns> /// <returns>resulting Polynomial</returns>
public static Polynomial operator -(Polynomial a, Polynomial b) public static Polynomial operator -(Polynomial a, Polynomial b)
{ {
return Substract(a, b); return Subtract(a, b);
} }
/// <summary> /// <summary>
@ -368,128 +365,147 @@ namespace MathNet.Numerics
} }
/// <summary> /// <summary>
/// pointwise division of two Polynomials /// Point-wise division of two Polynomials
/// </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 DividePointwise(Polynomial a, Polynomial b) public static Polynomial PointwiseDivide(Polynomial a, Polynomial b)
{ {
var aa = a.Clone() as Polynomial; var ac = a.Coefficients;
var bb = b.Clone() as Polynomial; 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; for (int i = commonLength; i < result.Length; i++)
double[] res = new double[aa.Coefficients.Length];
for (int ii = 0; ii < n; ii++)
{ {
res[ii] = aa.Coefficients[ii] / bb.Coefficients[ii]; result[i] = ac[i] / 0.0;
} }
return new Polynomial(res); return new Polynomial(result);
} }
/// <summary> /// <summary>
/// pointwise multiplication of two Polynomials /// Point-wise multiplication of two Polynomials
/// </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 MultiplyPointwise(Polynomial a, Polynomial b) public static Polynomial PointwiseMultiply(Polynomial a, Polynomial b)
{ {
var aa = a.Clone() as Polynomial; var ac = a.Coefficients;
var bb = b.Clone() as Polynomial; var bc = b.Coefficients;
if (aa.Coefficients.Length != bb.Coefficients.Length)
{
MakeSameLength(ref aa, ref bb);
}
int n = aa.Coefficients.Length; var degree = Math.Min(a.Degree, b.Degree);
double[] res = new double[aa.Coefficients.Length]; var result = new double[degree + 1];
for (int ii = 0; ii < n; ii++) 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);
} }
/// <summary> /// <summary>
/// Addition of two Polynomials (piecewise) /// Addition of two Polynomials (point-wise).
/// </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 Add(Polynomial a, Polynomial b) public static Polynomial Add(Polynomial a, Polynomial b)
{ {
var aa = a.Clone() as Polynomial; var ac = a.Coefficients;
var bb = b.Clone() as Polynomial; 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; int acLength = Math.Min(ac.Length, result.Length);
double[] res = new double[n]; for (int i = commonLength; i < acLength; i++)
for (int ii = 0; ii < n; ii++)
{ {
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);
} }
/// <summary> /// <summary>
/// substraction of two Polynomials (piecewise) /// Subtraction of two Polynomials (point-wise).
/// </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 Substract(Polynomial a, Polynomial b) public static Polynomial Subtract(Polynomial a, Polynomial b)
{ {
var aa = a.Clone() as Polynomial; var ac = a.Coefficients;
var bb = b.Clone() as Polynomial; 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; int bcLength = Math.Min(bc.Length, result.Length);
double[] res = new double[n]; for (int i = commonLength; i < bcLength; i++)
for (int ii = 0; ii < n; ii++)
{ {
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);
} }
/// <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="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> DivideLong(Polynomial a, Polynomial b)
{ {
if (a == null) if (a == null)
throw new ArgumentNullException("a"); throw new ArgumentNullException(nameof(a));
if (b == null) 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"); 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"); 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"); throw new DivideByZeroException("b polynomial ends with zero");
var c1 = a.Coefficients.ToArray(); var c1 = a.Coefficients.ToArray();
@ -572,8 +588,8 @@ namespace MathNet.Numerics
/// <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> DivideLong(Polynomial b)
{ {
return DivideLong(this, b); return DivideLong(this, b);
@ -643,30 +659,6 @@ namespace MathNet.Numerics
return Coefficients.ToArray(); 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);
}
}
/// <summary> /// <summary>
/// Full convolution of two arrays /// Full convolution of two arrays
/// </summary> /// </summary>

Write Preview
Loading…
Cancel
Save