Browse Source

Polynomial: migrate evaluation into Polynomial type

ridge-regression
Christoph Ruegg 8 years ago
parent
commit
9f8da8b00c
  1. 2
      src/FSharp/Fit.fs
  2. 2
      src/Numerics.Tests/FitTests.cs
  3. 4
      src/Numerics.Tests/RootFindingTests/BrentTest.cs
  4. 12
      src/Numerics.Tests/RootFindingTests/CubicTest.cs
  5. 8
      src/Numerics.Tests/RootFindingTests/NewtonRaphsonTest.cs
  6. 8
      src/Numerics.Tests/RootFindingTests/RobustNewtonRaphsonTest.cs
  7. 4
      src/Numerics.Tests/RootFindingTests/SecantTest.cs
  8. 8
      src/Numerics.Tests/SpecialFunctionsTests/ModifiedBesselTests.cs
  9. 6
      src/Numerics/Fit.cs
  10. 61
      src/Numerics/Polynomial.cs
  11. 42
      src/Numerics/SpecialFunctions/Erf.cs
  12. 30
      src/Numerics/SpecialFunctions/Evaluate.cs

2
src/FSharp/Fit.fs

@ -64,7 +64,7 @@ module Fit =
let multiDimFunc intercept x y = Fit.MultiDimFunc(x,y,intercept) |> tofs let multiDimFunc intercept x y = Fit.MultiDimFunc(x,y,intercept) |> tofs
/// Least-Squares fitting the points (x,y) to a k-order polynomial y : x -> p0 + p1*x + p2*x^2 + ... + pk*x^k, /// Least-Squares fitting the points (x,y) to a k-order polynomial y : x -> p0 + p1*x + p2*x^2 + ... + pk*x^k,
/// returning its best fitting parameters as [p0, p1, p2, ..., pk] array, compatible with Evaluate.Polynomial. /// returning its best fitting parameters as [p0, p1, p2, ..., pk] array, compatible with Polynomial.Evaluate.
let polynomial order x y = Fit.Polynomial(x,y,order) let polynomial order x y = Fit.Polynomial(x,y,order)
/// Least-Squares fitting the points (x,y) to a k-order polynomial y : x -> p0 + p1*x + p2*x^2 + ... + pk*x^k, /// Least-Squares fitting the points (x,y) to a k-order polynomial y : x -> p0 + p1*x + p2*x^2 + ... + pk*x^k,

2
src/Numerics.Tests/FitTests.cs

@ -273,7 +273,7 @@ namespace MathNet.Numerics.UnitTests
var resf = Fit.PolynomialFunc(x, y, 2); var resf = Fit.PolynomialFunc(x, y, 2);
foreach (var z in Enumerable.Range(-3, 10)) foreach (var z in Enumerable.Range(-3, 10))
{ {
Assert.AreEqual(Evaluate.Polynomial(z, resp), resf(z), 1e-4); Assert.AreEqual(Polynomial.Evaluate(z, resp), resf(z), 1e-4);
} }
} }

4
src/Numerics.Tests/RootFindingTests/BrentTest.cs

@ -70,11 +70,11 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests
public void Cubic() public void Cubic()
{ {
// with complex roots (looking for the real root only): 3x^3 + 4x^2 + 5x + 6 // with complex roots (looking for the real root only): 3x^3 + 4x^2 + 5x + 6
Func<double, double> f1 = x => Evaluate.Polynomial(x, 6, 5, 4, 3); Func<double, double> f1 = x => Polynomial.Evaluate(x, 6, 5, 4, 3);
Assert.AreEqual(-1.265328088928, Brent.FindRoot(f1, -2, -1, 1e-8, 100), 1e-6); Assert.AreEqual(-1.265328088928, Brent.FindRoot(f1, -2, -1, 1e-8, 100), 1e-6);
// real roots only: 2x^3 + 4x^2 - 50x + 6 // real roots only: 2x^3 + 4x^2 - 50x + 6
Func<double, double> f2 = x => Evaluate.Polynomial(x, 6, -50, 4, 2); Func<double, double> f2 = x => Polynomial.Evaluate(x, 6, -50, 4, 2);
Assert.AreEqual(-6.1466562197069, Brent.FindRoot(f2, -6.5, -5.5, 1e-8, 100), 1e-6); Assert.AreEqual(-6.1466562197069, Brent.FindRoot(f2, -6.5, -5.5, 1e-8, 100), 1e-6);
Assert.AreEqual(0.12124737195841, Brent.FindRoot(f2, -0.5, 0.5, 1e-8, 100), 1e-6); Assert.AreEqual(0.12124737195841, Brent.FindRoot(f2, -0.5, 0.5, 1e-8, 100), 1e-6);
Assert.AreEqual(4.0254088477485, Brent.FindRoot(f2, 3.5, 4.5, 1e-8, 100), 1e-6); Assert.AreEqual(4.0254088477485, Brent.FindRoot(f2, 3.5, 4.5, 1e-8, 100), 1e-6);

12
src/Numerics.Tests/RootFindingTests/CubicTest.cs

@ -106,12 +106,12 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests
public void ComplexRootsAreRoots(double d, double c, double b, double a) public void ComplexRootsAreRoots(double d, double c, double b, double a)
{ {
var roots = FindRoots.Cubic(d, c, b, a); var roots = FindRoots.Cubic(d, c, b, a);
Assert.That(Evaluate.Polynomial(roots.Item1, d, c, b, a).Real, Is.EqualTo(0).Within(1e-12)); Assert.That(Polynomial.Evaluate(roots.Item1, d, c, b, a).Real, Is.EqualTo(0).Within(1e-12));
Assert.That(Evaluate.Polynomial(roots.Item1, d, c, b, a).Imaginary, Is.EqualTo(0).Within(1e-12)); Assert.That(Polynomial.Evaluate(roots.Item1, d, c, b, a).Imaginary, Is.EqualTo(0).Within(1e-12));
Assert.That(Evaluate.Polynomial(roots.Item2, d, c, b, a).Real, Is.EqualTo(0).Within(1e-12)); Assert.That(Polynomial.Evaluate(roots.Item2, d, c, b, a).Real, Is.EqualTo(0).Within(1e-12));
Assert.That(Evaluate.Polynomial(roots.Item2, d, c, b, a).Imaginary, Is.EqualTo(0).Within(1e-12)); Assert.That(Polynomial.Evaluate(roots.Item2, d, c, b, a).Imaginary, Is.EqualTo(0).Within(1e-12));
Assert.That(Evaluate.Polynomial(roots.Item3, d, c, b, a).Real, Is.EqualTo(0).Within(1e-12)); Assert.That(Polynomial.Evaluate(roots.Item3, d, c, b, a).Real, Is.EqualTo(0).Within(1e-12));
Assert.That(Evaluate.Polynomial(roots.Item3, d, c, b, a).Imaginary, Is.EqualTo(0).Within(1e-12)); Assert.That(Polynomial.Evaluate(roots.Item3, d, c, b, a).Imaginary, Is.EqualTo(0).Within(1e-12));
} }
} }
} }

8
src/Numerics.Tests/RootFindingTests/NewtonRaphsonTest.cs

@ -95,14 +95,14 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests
public void Cubic() public void Cubic()
{ {
// with complex roots (looking for the real root only): 3x^3 + 4x^2 + 5x + 6, derivative 9x^2 + 8x + 5 // with complex roots (looking for the real root only): 3x^3 + 4x^2 + 5x + 6, derivative 9x^2 + 8x + 5
Func<double, double> f1 = x => Evaluate.Polynomial(x, 6, 5, 4, 3); Func<double, double> f1 = x => Polynomial.Evaluate(x, 6, 5, 4, 3);
Func<double, double> df1 = x => Evaluate.Polynomial(x, 5, 8, 9); Func<double, double> df1 = x => Polynomial.Evaluate(x, 5, 8, 9);
Assert.AreEqual(-1.265328088928, NewtonRaphson.FindRoot(f1, df1, -2, -1, 1e-10), 1e-6); Assert.AreEqual(-1.265328088928, NewtonRaphson.FindRoot(f1, df1, -2, -1, 1e-10), 1e-6);
Assert.AreEqual(-1.265328088928, NewtonRaphson.FindRoot(f1, df1, -5, 5, 1e-10), 1e-6); Assert.AreEqual(-1.265328088928, NewtonRaphson.FindRoot(f1, df1, -5, 5, 1e-10), 1e-6);
// real roots only: 2x^3 + 4x^2 - 50x + 6, derivative 6x^2 + 8x - 50 // real roots only: 2x^3 + 4x^2 - 50x + 6, derivative 6x^2 + 8x - 50
Func<double, double> f2 = x => Evaluate.Polynomial(x, 6, -50, 4, 2); Func<double, double> f2 = x => Polynomial.Evaluate(x, 6, -50, 4, 2);
Func<double, double> df2 = x => Evaluate.Polynomial(x, -50, 8, 6); Func<double, double> df2 = x => Polynomial.Evaluate(x, -50, 8, 6);
Assert.AreEqual(-6.1466562197069, NewtonRaphson.FindRoot(f2, df2, -8, -5, 1e-10), 1e-6); Assert.AreEqual(-6.1466562197069, NewtonRaphson.FindRoot(f2, df2, -8, -5, 1e-10), 1e-6);
Assert.AreEqual(0.12124737195841, NewtonRaphson.FindRoot(f2, df2, -1, 1, 1e-10), 1e-6); Assert.AreEqual(0.12124737195841, NewtonRaphson.FindRoot(f2, df2, -1, 1, 1e-10), 1e-6);
Assert.AreEqual(4.0254088477485, NewtonRaphson.FindRoot(f2, df2, 3, 5, 1e-10), 1e-6); Assert.AreEqual(4.0254088477485, NewtonRaphson.FindRoot(f2, df2, 3, 5, 1e-10), 1e-6);

8
src/Numerics.Tests/RootFindingTests/RobustNewtonRaphsonTest.cs

@ -102,14 +102,14 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests
public void Cubic() public void Cubic()
{ {
// with complex roots (looking for the real root only): 3x^3 + 4x^2 + 5x + 6, derivative 9x^2 + 8x + 5 // with complex roots (looking for the real root only): 3x^3 + 4x^2 + 5x + 6, derivative 9x^2 + 8x + 5
Func<double, double> f1 = x => Evaluate.Polynomial(x, 6, 5, 4, 3); Func<double, double> f1 = x => Polynomial.Evaluate(x, 6, 5, 4, 3);
Func<double, double> df1 = x => Evaluate.Polynomial(x, 5, 8, 9); Func<double, double> df1 = x => Polynomial.Evaluate(x, 5, 8, 9);
Assert.AreEqual(-1.265328088928, RobustNewtonRaphson.FindRoot(f1, df1, -2, -1, 1e-10, 100, 20), 1e-6); Assert.AreEqual(-1.265328088928, RobustNewtonRaphson.FindRoot(f1, df1, -2, -1, 1e-10, 100, 20), 1e-6);
Assert.AreEqual(-1.265328088928, RobustNewtonRaphson.FindRoot(f1, df1, -5, 5, 1e-10, 100, 20), 1e-6); Assert.AreEqual(-1.265328088928, RobustNewtonRaphson.FindRoot(f1, df1, -5, 5, 1e-10, 100, 20), 1e-6);
// real roots only: 2x^3 + 4x^2 - 50x + 6, derivative 6x^2 + 8x - 50 // real roots only: 2x^3 + 4x^2 - 50x + 6, derivative 6x^2 + 8x - 50
Func<double, double> f2 = x => Evaluate.Polynomial(x, 6, -50, 4, 2); Func<double, double> f2 = x => Polynomial.Evaluate(x, 6, -50, 4, 2);
Func<double, double> df2 = x => Evaluate.Polynomial(x, -50, 8, 6); Func<double, double> df2 = x => Polynomial.Evaluate(x, -50, 8, 6);
Assert.AreEqual(-6.1466562197069, RobustNewtonRaphson.FindRoot(f2, df2, -8, -5, 1e-10, 100, 20), 1e-6); Assert.AreEqual(-6.1466562197069, RobustNewtonRaphson.FindRoot(f2, df2, -8, -5, 1e-10, 100, 20), 1e-6);
Assert.AreEqual(0.12124737195841, RobustNewtonRaphson.FindRoot(f2, df2, -1, 1, 1e-10, 100, 20), 1e-6); Assert.AreEqual(0.12124737195841, RobustNewtonRaphson.FindRoot(f2, df2, -1, 1, 1e-10, 100, 20), 1e-6);
Assert.AreEqual(4.0254088477485, RobustNewtonRaphson.FindRoot(f2, df2, 3, 5, 1e-10, 100, 20), 1e-6); Assert.AreEqual(4.0254088477485, RobustNewtonRaphson.FindRoot(f2, df2, 3, 5, 1e-10, 100, 20), 1e-6);

4
src/Numerics.Tests/RootFindingTests/SecantTest.cs

@ -68,12 +68,12 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests
public void Cubic() public void Cubic()
{ {
// with complex roots (looking for the real root only): 3x^3 + 4x^2 + 5x + 6 // with complex roots (looking for the real root only): 3x^3 + 4x^2 + 5x + 6
Func<double, double> f1 = x => Evaluate.Polynomial(x, 6, 5, 4, 3); Func<double, double> f1 = x => Polynomial.Evaluate(x, 6, 5, 4, 3);
Assert.AreEqual(-1.265328088928, Secant.FindRoot(f1, -2, -1, -10, 10, 1e-10), 1e-6); Assert.AreEqual(-1.265328088928, Secant.FindRoot(f1, -2, -1, -10, 10, 1e-10), 1e-6);
Assert.AreEqual(-1.265328088928, Secant.FindRoot(f1, -5, 6, -10, 10, 1e-10), 1e-6); Assert.AreEqual(-1.265328088928, Secant.FindRoot(f1, -5, 6, -10, 10, 1e-10), 1e-6);
// real roots only: 2x^3 + 4x^2 - 50x + 6 // real roots only: 2x^3 + 4x^2 - 50x + 6
Func<double, double> f2 = x => Evaluate.Polynomial(x, 6, -50, 4, 2); Func<double, double> f2 = x => Polynomial.Evaluate(x, 6, -50, 4, 2);
Assert.AreEqual(-6.1466562197069, Secant.FindRoot(f2, -8, -5, -10, 10, 1e-10), 1e-6); Assert.AreEqual(-6.1466562197069, Secant.FindRoot(f2, -8, -5, -10, 10, 1e-10), 1e-6);
Assert.AreEqual(0.12124737195841, Secant.FindRoot(f2, -1, 1, -10, 10, 1e-10), 1e-6); Assert.AreEqual(0.12124737195841, Secant.FindRoot(f2, -1, 1, -10, 10, 1e-10), 1e-6);
Assert.AreEqual(4.0254088477485, Secant.FindRoot(f2, 3, 5, 0, 10, 1e-10), 1e-6); Assert.AreEqual(4.0254088477485, Secant.FindRoot(f2, 3, 5, 0, 10, 1e-10), 1e-6);

8
src/Numerics.Tests/SpecialFunctionsTests/ModifiedBesselTests.cs

@ -42,7 +42,7 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests
public void BesselI0Approx([Range(-3.75, 3.75, 0.25)] double x) public void BesselI0Approx([Range(-3.75, 3.75, 0.25)] double x)
{ {
// Approx by Abramowitz/Stegun 9.8.1 // Approx by Abramowitz/Stegun 9.8.1
Assert.AreEqual(Evaluate.Polynomial(x/3.75, 1.0, 0.0, 3.5156229, 0.0, 3.0899424, 0.0, 1.2067492, 0.0, 0.2659732, 0.0, 0.0360768, 0.0, 0.0045813), SpecialFunctions.BesselI0(x), 1e-7); Assert.AreEqual(Polynomial.Evaluate(x/3.75, 1.0, 0.0, 3.5156229, 0.0, 3.0899424, 0.0, 1.2067492, 0.0, 0.2659732, 0.0, 0.0360768, 0.0, 0.0045813), SpecialFunctions.BesselI0(x), 1e-7);
} }
[TestCase(0.0, 1.0)] [TestCase(0.0, 1.0)]
@ -62,7 +62,7 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests
public void BesselI1Approx([Range(-3.75, 3.75, 0.25)] double x) public void BesselI1Approx([Range(-3.75, 3.75, 0.25)] double x)
{ {
// Approx by Abramowitz/Stegun 9.8.3 // Approx by Abramowitz/Stegun 9.8.3
Assert.AreEqual(Evaluate.Polynomial(x/3.75, 0.5, 0.0, 0.87890594, 0.0, 0.51498869, 0.0, 0.15084934, 0.0, 0.02658733, 0.0, 0.00301532, 0.0, 0.00032411)*x, SpecialFunctions.BesselI1(x), 1e-8); Assert.AreEqual(Polynomial.Evaluate(x/3.75, 0.5, 0.0, 0.87890594, 0.0, 0.51498869, 0.0, 0.15084934, 0.0, 0.02658733, 0.0, 0.00301532, 0.0, 0.00032411)*x, SpecialFunctions.BesselI1(x), 1e-8);
} }
[TestCase(0.0, 0.0)] [TestCase(0.0, 0.0)]
@ -82,7 +82,7 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests
public void BesselK0Approx([Range(0.20, 2.0, 0.20)] double x) public void BesselK0Approx([Range(0.20, 2.0, 0.20)] double x)
{ {
// Approx by Abramowitz/Stegun 9.8.5 // Approx by Abramowitz/Stegun 9.8.5
Assert.AreEqual(Evaluate.Polynomial(x/2.0, -Math.Log(x/2.0)*SpecialFunctions.BesselI0(x) - 0.57721566, 0.0, 0.42278420, 0.0, 0.23069756, 0.0, 0.03488590, 0.0, 0.00262698, 0.0, 0.00010750, 0.0, 0.00000740), SpecialFunctions.BesselK0(x), 1e-8); Assert.AreEqual(Polynomial.Evaluate(x/2.0, -Math.Log(x/2.0)*SpecialFunctions.BesselI0(x) - 0.57721566, 0.0, 0.42278420, 0.0, 0.23069756, 0.0, 0.03488590, 0.0, 0.00262698, 0.0, 0.00010750, 0.0, 0.00000740), SpecialFunctions.BesselK0(x), 1e-8);
} }
[TestCase(1e-10, 23.14178244559887)] [TestCase(1e-10, 23.14178244559887)]
@ -101,7 +101,7 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests
public void BesselK1Approx([Range(0.20, 2.0, 0.20)] double x) public void BesselK1Approx([Range(0.20, 2.0, 0.20)] double x)
{ {
// Approx by Abramowitz/Stegun 9.8.7 // Approx by Abramowitz/Stegun 9.8.7
Assert.AreEqual(Evaluate.Polynomial(x/2.0, x*Math.Log(x/2.0)*SpecialFunctions.BesselI1(x) + 1.0, 0.0, 0.15443144, 0.0, -0.67278579, 0.0, -0.18156897, 0.0, -0.01919402, 0.0, -0.00110404, 0.0, -0.00004686), SpecialFunctions.BesselK1(x)*x, 1e-8); Assert.AreEqual(Polynomial.Evaluate(x/2.0, x*Math.Log(x/2.0)*SpecialFunctions.BesselI1(x) + 1.0, 0.0, 0.15443144, 0.0, -0.67278579, 0.0, -0.18156897, 0.0, -0.01919402, 0.0, -0.00110404, 0.0, -0.00004686), SpecialFunctions.BesselK1(x)*x, 1e-8);
} }
[TestCase(1e-10, 1.0e+10)] [TestCase(1e-10, 1.0e+10)]

6
src/Numerics/Fit.cs

@ -154,7 +154,7 @@ namespace MathNet.Numerics
/// <summary> /// <summary>
/// Least-Squares fitting the points (x,y) to a k-order polynomial y : x -> p0 + p1*x + p2*x^2 + ... + pk*x^k, /// Least-Squares fitting the points (x,y) to a k-order polynomial y : x -> p0 + p1*x + p2*x^2 + ... + pk*x^k,
/// returning its best fitting parameters as [p0, p1, p2, ..., pk] array, compatible with Evaluate.Polynomial. /// returning its best fitting parameters as [p0, p1, p2, ..., pk] array, compatible with Polynomial.Evaluate.
/// A polynomial with order/degree k has (k+1) coefficients and thus requires at least (k+1) samples. /// A polynomial with order/degree k has (k+1) coefficients and thus requires at least (k+1) samples.
/// </summary> /// </summary>
public static double[] Polynomial(double[] x, double[] y, int order, DirectRegressionMethod method = DirectRegressionMethod.QR) public static double[] Polynomial(double[] x, double[] y, int order, DirectRegressionMethod method = DirectRegressionMethod.QR)
@ -171,12 +171,12 @@ namespace MathNet.Numerics
public static Func<double, double> PolynomialFunc(double[] x, double[] y, int order, DirectRegressionMethod method = DirectRegressionMethod.QR) public static Func<double, double> PolynomialFunc(double[] x, double[] y, int order, DirectRegressionMethod method = DirectRegressionMethod.QR)
{ {
var parameters = Polynomial(x, y, order, method); var parameters = Polynomial(x, y, order, method);
return z => Evaluate.Polynomial(z, parameters); return z => Numerics.Polynomial.Evaluate(z, parameters);
} }
/// <summary> /// <summary>
/// Weighted Least-Squares fitting the points (x,y) and weights w to a k-order polynomial y : x -> p0 + p1*x + p2*x^2 + ... + pk*x^k, /// Weighted Least-Squares fitting the points (x,y) and weights w to a k-order polynomial y : x -> p0 + p1*x + p2*x^2 + ... + pk*x^k,
/// returning its best fitting parameters as [p0, p1, p2, ..., pk] array, compatible with Evaluate.Polynomial. /// returning its best fitting parameters as [p0, p1, p2, ..., pk] array, compatible with Polynomial.Evaluate.
/// A polynomial with order/degree k has (k+1) coefficients and thus requires at least (k+1) samples. /// A polynomial with order/degree k has (k+1) coefficients and thus requires at least (k+1) samples.
/// </summary> /// </summary>
public static double[] PolynomialWeighted(double[] x, double[] y, double[] w, int order) public static double[] PolynomialWeighted(double[] x, double[] y, double[] w, int order)

61
src/Numerics/Polynomial.cs

@ -119,13 +119,70 @@ namespace MathNet.Numerics
} }
#region Evaluation #region Evaluation
/// <summary>
/// Evaluate a polynomial at point x.
/// Coefficients are ordered by power with power k at index k.
/// Example: coefficients [3,-1,2] represent y=2x^2-x+3.
/// </summary>
/// <param name="z">The location where to evaluate the polynomial at.</param>
/// <param name="coefficients">The coefficients of the polynomial, coefficient for power k at index k.</param>
public static double Evaluate(double z, params double[] coefficients)
{
double sum = coefficients[coefficients.Length - 1];
for (int i = coefficients.Length - 2; i >= 0; --i)
{
sum *= z;
sum += coefficients[i];
}
return sum;
}
/// <summary>
/// Evaluate a polynomial at point x.
/// Coefficients are ordered by power with power k at index k.
/// Example: coefficients [3,-1,2] represent y=2x^2-x+3.
/// </summary>
/// <param name="z">The location where to evaluate the polynomial at.</param>
/// <param name="coefficients">The coefficients of the polynomial, coefficient for power k at index k.</param>
public static Complex Evaluate(Complex z, params double[] coefficients)
{
Complex sum = coefficients[coefficients.Length - 1];
for (int i = coefficients.Length - 2; i >= 0; --i)
{
sum *= z;
sum += coefficients[i];
}
return sum;
}
/// <summary>
/// Evaluate a polynomial at point x.
/// Coefficients are ordered by power with power k at index k.
/// Example: coefficients [3,-1,2] represent y=2x^2-x+3.
/// </summary>
/// <param name="z">The location where to evaluate the polynomial at.</param>
/// <param name="coefficients">The coefficients of the polynomial, coefficient for power k at index k.</param>
public static Complex Evaluate(Complex z, params Complex[] coefficients)
{
Complex sum = coefficients[coefficients.Length - 1];
for (int i = coefficients.Length - 2; i >= 0; --i)
{
sum *= z;
sum += coefficients[i];
}
return sum;
}
/// <summary> /// <summary>
/// Evaluate a polynomial at point x. /// Evaluate a polynomial at point x.
/// </summary> /// </summary>
/// <param name="z">The location where to evaluate the polynomial at.</param> /// <param name="z">The location where to evaluate the polynomial at.</param>
public double Evaluate(double z) public double Evaluate(double z)
{ {
return Numerics.Evaluate.Polynomial(z, Coefficients); return Evaluate(z, Coefficients);
} }
/// <summary> /// <summary>
@ -134,7 +191,7 @@ namespace MathNet.Numerics
/// <param name="z">The location where to evaluate the polynomial at.</param> /// <param name="z">The location where to evaluate the polynomial at.</param>
public Complex Evaluate(Complex z) public Complex Evaluate(Complex z)
{ {
return Numerics.Evaluate.Polynomial(z, Coefficients); return Evaluate(z, Coefficients);
} }
/// <summary> /// <summary>

42
src/Numerics/SpecialFunctions/Erf.cs

@ -418,7 +418,7 @@ namespace MathNet.Numerics
else else
{ {
// Worst case absolute error found: 6.688618532e-21 // Worst case absolute error found: 6.688618532e-21
result = (z*1.125) + (z*Evaluate.Polynomial(z, ErfImpAn)/Evaluate.Polynomial(z, ErfImpAd)); result = (z*1.125) + (z*Polynomial.Evaluate(z, ErfImpAn)/Polynomial.Evaluate(z, ErfImpAd));
} }
} }
else if (z < 110) else if (z < 110)
@ -429,79 +429,79 @@ namespace MathNet.Numerics
if (z < 0.75) if (z < 0.75)
{ {
// Worst case absolute error found: 5.582813374e-21 // Worst case absolute error found: 5.582813374e-21
r = Evaluate.Polynomial(z - 0.5, ErfImpBn)/Evaluate.Polynomial(z - 0.5, ErfImpBd); r = Polynomial.Evaluate(z - 0.5, ErfImpBn)/Polynomial.Evaluate(z - 0.5, ErfImpBd);
b = 0.3440242112F; b = 0.3440242112F;
} }
else if (z < 1.25) else if (z < 1.25)
{ {
// Worst case absolute error found: 4.01854729e-21 // Worst case absolute error found: 4.01854729e-21
r = Evaluate.Polynomial(z - 0.75, ErfImpCn)/Evaluate.Polynomial(z - 0.75, ErfImpCd); r = Polynomial.Evaluate(z - 0.75, ErfImpCn)/Polynomial.Evaluate(z - 0.75, ErfImpCd);
b = 0.419990927F; b = 0.419990927F;
} }
else if (z < 2.25) else if (z < 2.25)
{ {
// Worst case absolute error found: 2.866005373e-21 // Worst case absolute error found: 2.866005373e-21
r = Evaluate.Polynomial(z - 1.25, ErfImpDn)/Evaluate.Polynomial(z - 1.25, ErfImpDd); r = Polynomial.Evaluate(z - 1.25, ErfImpDn)/Polynomial.Evaluate(z - 1.25, ErfImpDd);
b = 0.4898625016F; b = 0.4898625016F;
} }
else if (z < 3.5) else if (z < 3.5)
{ {
// Worst case absolute error found: 1.045355789e-21 // Worst case absolute error found: 1.045355789e-21
r = Evaluate.Polynomial(z - 2.25, ErfImpEn) /Evaluate.Polynomial(z - 2.25, ErfImpEd); r = Polynomial.Evaluate(z - 2.25, ErfImpEn) /Polynomial.Evaluate(z - 2.25, ErfImpEd);
b = 0.5317370892F; b = 0.5317370892F;
} }
else if (z < 5.25) else if (z < 5.25)
{ {
// Worst case absolute error found: 8.300028706e-22 // Worst case absolute error found: 8.300028706e-22
r = Evaluate.Polynomial(z - 3.5, ErfImpFn) /Evaluate.Polynomial(z - 3.5, ErfImpFd); r = Polynomial.Evaluate(z - 3.5, ErfImpFn) /Polynomial.Evaluate(z - 3.5, ErfImpFd);
b = 0.5489973426F; b = 0.5489973426F;
} }
else if (z < 8) else if (z < 8)
{ {
// Worst case absolute error found: 1.700157534e-21 // Worst case absolute error found: 1.700157534e-21
r = Evaluate.Polynomial(z - 5.25, ErfImpGn) /Evaluate.Polynomial(z - 5.25, ErfImpGd); r = Polynomial.Evaluate(z - 5.25, ErfImpGn) /Polynomial.Evaluate(z - 5.25, ErfImpGd);
b = 0.5571740866F; b = 0.5571740866F;
} }
else if (z < 11.5) else if (z < 11.5)
{ {
// Worst case absolute error found: 3.002278011e-22 // Worst case absolute error found: 3.002278011e-22
r = Evaluate.Polynomial(z - 8, ErfImpHn) /Evaluate.Polynomial(z - 8, ErfImpHd); r = Polynomial.Evaluate(z - 8, ErfImpHn) /Polynomial.Evaluate(z - 8, ErfImpHd);
b = 0.5609807968F; b = 0.5609807968F;
} }
else if (z < 17) else if (z < 17)
{ {
// Worst case absolute error found: 6.741114695e-21 // Worst case absolute error found: 6.741114695e-21
r = Evaluate.Polynomial(z - 11.5, ErfImpIn) /Evaluate.Polynomial(z - 11.5, ErfImpId); r = Polynomial.Evaluate(z - 11.5, ErfImpIn) /Polynomial.Evaluate(z - 11.5, ErfImpId);
b = 0.5626493692F; b = 0.5626493692F;
} }
else if (z < 24) else if (z < 24)
{ {
// Worst case absolute error found: 7.802346984e-22 // Worst case absolute error found: 7.802346984e-22
r = Evaluate.Polynomial(z - 17, ErfImpJn) /Evaluate.Polynomial(z - 17, ErfImpJd); r = Polynomial.Evaluate(z - 17, ErfImpJn) /Polynomial.Evaluate(z - 17, ErfImpJd);
b = 0.5634598136F; b = 0.5634598136F;
} }
else if (z < 38) else if (z < 38)
{ {
// Worst case absolute error found: 2.414228989e-22 // Worst case absolute error found: 2.414228989e-22
r = Evaluate.Polynomial(z - 24, ErfImpKn) /Evaluate.Polynomial(z - 24, ErfImpKd); r = Polynomial.Evaluate(z - 24, ErfImpKn) /Polynomial.Evaluate(z - 24, ErfImpKd);
b = 0.5638477802F; b = 0.5638477802F;
} }
else if (z < 60) else if (z < 60)
{ {
// Worst case absolute error found: 5.896543869e-24 // Worst case absolute error found: 5.896543869e-24
r = Evaluate.Polynomial(z - 38, ErfImpLn) /Evaluate.Polynomial(z - 38, ErfImpLd); r = Polynomial.Evaluate(z - 38, ErfImpLn) /Polynomial.Evaluate(z - 38, ErfImpLd);
b = 0.5640528202F; b = 0.5640528202F;
} }
else if (z < 85) else if (z < 85)
{ {
// Worst case absolute error found: 3.080612264e-21 // Worst case absolute error found: 3.080612264e-21
r = Evaluate.Polynomial(z - 60, ErfImpMn) /Evaluate.Polynomial(z - 60, ErfImpMd); r = Polynomial.Evaluate(z - 60, ErfImpMn) /Polynomial.Evaluate(z - 60, ErfImpMd);
b = 0.5641309023F; b = 0.5641309023F;
} }
else else
{ {
// Worst case absolute error found: 8.094633491e-22 // Worst case absolute error found: 8.094633491e-22
r = Evaluate.Polynomial(z - 85, ErfImpNn) /Evaluate.Polynomial(z - 85, ErfImpNd); r = Polynomial.Evaluate(z - 85, ErfImpNn) /Polynomial.Evaluate(z - 85, ErfImpNd);
b = 0.5641584396F; b = 0.5641584396F;
} }
@ -589,7 +589,7 @@ namespace MathNet.Numerics
// Maximum Deviation Found (actual error term at infinite precision) 8.030e-21 // Maximum Deviation Found (actual error term at infinite precision) 8.030e-21
const float y = 0.0891314744949340820313f; const float y = 0.0891314744949340820313f;
double g = p*(p + 10); double g = p*(p + 10);
double r = Evaluate.Polynomial(p, ErvInvImpAn) /Evaluate.Polynomial(p, ErvInvImpAd); double r = Polynomial.Evaluate(p, ErvInvImpAn) /Polynomial.Evaluate(p, ErvInvImpAd);
result = (g*y) + (g*r); result = (g*y) + (g*r);
} }
else if (q >= 0.25) else if (q >= 0.25)
@ -607,7 +607,7 @@ namespace MathNet.Numerics
const float y = 2.249481201171875f; const float y = 2.249481201171875f;
double g = Math.Sqrt(-2*Math.Log(q)); double g = Math.Sqrt(-2*Math.Log(q));
double xs = q - 0.25; double xs = q - 0.25;
double r = Evaluate.Polynomial(xs, ErvInvImpBn) /Evaluate.Polynomial(xs, ErvInvImpBd); double r = Polynomial.Evaluate(xs, ErvInvImpBn) /Polynomial.Evaluate(xs, ErvInvImpBd);
result = g/(y + r); result = g/(y + r);
} }
else else
@ -635,7 +635,7 @@ namespace MathNet.Numerics
// Max error found: 1.089051e-20 // Max error found: 1.089051e-20
const float y = 0.807220458984375f; const float y = 0.807220458984375f;
double xs = x - 1.125; double xs = x - 1.125;
double r = Evaluate.Polynomial(xs, ErvInvImpCn) /Evaluate.Polynomial(xs, ErvInvImpCd); double r = Polynomial.Evaluate(xs, ErvInvImpCn) /Polynomial.Evaluate(xs, ErvInvImpCd);
result = (y*x) + (r*x); result = (y*x) + (r*x);
} }
else if (x < 6) else if (x < 6)
@ -643,7 +643,7 @@ namespace MathNet.Numerics
// Max error found: 8.389174e-21 // Max error found: 8.389174e-21
const float y = 0.93995571136474609375f; const float y = 0.93995571136474609375f;
double xs = x - 3; double xs = x - 3;
double r = Evaluate.Polynomial(xs, ErvInvImpDn) /Evaluate.Polynomial(xs, ErvInvImpDd); double r = Polynomial.Evaluate(xs, ErvInvImpDn) /Polynomial.Evaluate(xs, ErvInvImpDd);
result = (y*x) + (r*x); result = (y*x) + (r*x);
} }
else if (x < 18) else if (x < 18)
@ -651,7 +651,7 @@ namespace MathNet.Numerics
// Max error found: 1.481312e-19 // Max error found: 1.481312e-19
const float y = 0.98362827301025390625f; const float y = 0.98362827301025390625f;
double xs = x - 6; double xs = x - 6;
double r = Evaluate.Polynomial(xs, ErvInvImpEn) /Evaluate.Polynomial(xs, ErvInvImpEd); double r = Polynomial.Evaluate(xs, ErvInvImpEn) /Polynomial.Evaluate(xs, ErvInvImpEd);
result = (y*x) + (r*x); result = (y*x) + (r*x);
} }
else if (x < 44) else if (x < 44)
@ -659,7 +659,7 @@ namespace MathNet.Numerics
// Max error found: 5.697761e-20 // Max error found: 5.697761e-20
const float y = 0.99714565277099609375f; const float y = 0.99714565277099609375f;
double xs = x - 18; double xs = x - 18;
double r = Evaluate.Polynomial(xs, ErvInvImpFn) /Evaluate.Polynomial(xs, ErvInvImpFd); double r = Polynomial.Evaluate(xs, ErvInvImpFn) /Polynomial.Evaluate(xs, ErvInvImpFd);
result = (y*x) + (r*x); result = (y*x) + (r*x);
} }
else else
@ -667,7 +667,7 @@ namespace MathNet.Numerics
// Max error found: 1.279746e-20 // Max error found: 1.279746e-20
const float y = 0.99941349029541015625f; const float y = 0.99941349029541015625f;
double xs = x - 44; double xs = x - 44;
double r = Evaluate.Polynomial(xs, ErvInvImpGn) /Evaluate.Polynomial(xs, ErvInvImpGd); double r = Polynomial.Evaluate(xs, ErvInvImpGn) /Polynomial.Evaluate(xs, ErvInvImpGd);
result = (y*x) + (r*x); result = (y*x) + (r*x);
} }
} }

30
src/Numerics/SpecialFunctions/Evaluate.cs

@ -62,16 +62,10 @@ namespace MathNet.Numerics
/// </summary> /// </summary>
/// <param name="z">The location where to evaluate the polynomial at.</param> /// <param name="z">The location where to evaluate the polynomial at.</param>
/// <param name="coefficients">The coefficients of the polynomial, coefficient for power k at index k.</param> /// <param name="coefficients">The coefficients of the polynomial, coefficient for power k at index k.</param>
[Obsolete("Use Polynomial.Evaluate instead. Will be removed in the next major version.")]
public static double Polynomial(double z, params double[] coefficients) public static double Polynomial(double z, params double[] coefficients)
{ {
double sum = coefficients[coefficients.Length - 1]; return Numerics.Polynomial.Evaluate(z, coefficients);
for (int i = coefficients.Length - 2; i >= 0; --i)
{
sum *= z;
sum += coefficients[i];
}
return sum;
} }
/// <summary> /// <summary>
@ -81,16 +75,10 @@ namespace MathNet.Numerics
/// </summary> /// </summary>
/// <param name="z">The location where to evaluate the polynomial at.</param> /// <param name="z">The location where to evaluate the polynomial at.</param>
/// <param name="coefficients">The coefficients of the polynomial, coefficient for power k at index k.</param> /// <param name="coefficients">The coefficients of the polynomial, coefficient for power k at index k.</param>
[Obsolete("Use Polynomial.Evaluate instead. Will be removed in the next major version.")]
public static Complex Polynomial(Complex z, params double[] coefficients) public static Complex Polynomial(Complex z, params double[] coefficients)
{ {
Complex sum = coefficients[coefficients.Length - 1]; return Numerics.Polynomial.Evaluate(z, coefficients);
for (int i = coefficients.Length - 2; i >= 0; --i)
{
sum *= z;
sum += coefficients[i];
}
return sum;
} }
/// <summary> /// <summary>
@ -100,16 +88,10 @@ namespace MathNet.Numerics
/// </summary> /// </summary>
/// <param name="z">The location where to evaluate the polynomial at.</param> /// <param name="z">The location where to evaluate the polynomial at.</param>
/// <param name="coefficients">The coefficients of the polynomial, coefficient for power k at index k.</param> /// <param name="coefficients">The coefficients of the polynomial, coefficient for power k at index k.</param>
[Obsolete("Use Polynomial.Evaluate instead. Will be removed in the next major version.")]
public static Complex Polynomial(Complex z, params Complex[] coefficients) public static Complex Polynomial(Complex z, params Complex[] coefficients)
{ {
Complex sum = coefficients[coefficients.Length - 1]; return Numerics.Polynomial.Evaluate(z, coefficients);
for (int i = coefficients.Length - 2; i >= 0; --i)
{
sum *= z;
sum += coefficients[i];
}
return sum;
} }
/// <summary> /// <summary>

Loading…
Cancel
Save