// <copyright file="Evaluate.cs" company="Math.NET">
|
// Math.NET Numerics, part of the Math.NET Project
|
// http://numerics.mathdotnet.com
|
// http://github.com/mathnet/mathnet-numerics
|
//
|
// Copyright (c) 2009-2013 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.
|
// </copyright>
|
|
// <contribution>
|
// CERN - European Laboratory for Particle Physics
|
// http://www.docjar.com/html/api/cern/jet/math/Bessel.java.html
|
// Copyright 1999 CERN - European Laboratory for Particle Physics.
|
// Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
|
// is hereby granted without fee, provided that the above copyright notice appear in all copies and
|
// that both that copyright notice and this permission notice appear in supporting documentation.
|
// CERN makes no representations about the suitability of this software for any purpose.
|
// It is provided "as is" without expressed or implied warranty.
|
// TOMS757 - Uncommon Special Functions (Fortran77) by Allan McLeod
|
// http://people.sc.fsu.edu/~jburkardt/f77_src/toms757/toms757.html
|
// Wei Wu
|
// Cephes Math Library, Stephen L. Moshier
|
// ALGLIB 2.0.1, Sergey Bochkanov
|
// </contribution>
|
|
using System;
|
using Complex = System.Numerics.Complex;
|
|
// ReSharper disable CheckNamespace
|
namespace IStation.Numerics
|
// ReSharper restore CheckNamespace
|
{
|
/// <summary>
|
/// Evaluation functions, useful for function approximation.
|
/// </summary>
|
public static class Evaluate
|
{
|
/// <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>
|
[Obsolete("Use Polynomial.Evaluate instead. Will be removed in the next major version.")]
|
public static double Polynomial(double z, params double[] coefficients)
|
{
|
return Numerics.Polynomial.Evaluate(z, coefficients);
|
}
|
|
/// <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>
|
[Obsolete("Use Polynomial.Evaluate instead. Will be removed in the next major version.")]
|
public static Complex Polynomial(Complex z, params double[] coefficients)
|
{
|
return Numerics.Polynomial.Evaluate(z, coefficients);
|
}
|
|
/// <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>
|
[Obsolete("Use Polynomial.Evaluate instead. Will be removed in the next major version.")]
|
public static Complex Polynomial(Complex z, params Complex[] coefficients)
|
{
|
return Numerics.Polynomial.Evaluate(z, coefficients);
|
}
|
|
/// <summary>
|
/// Numerically stable series summation
|
/// </summary>
|
/// <param name="nextSummand">provides the summands sequentially</param>
|
/// <returns>Sum</returns>
|
internal static double Series(Func<double> nextSummand)
|
{
|
double compensation = 0.0;
|
double current;
|
const double factor = 1 << 16;
|
|
double sum = nextSummand();
|
|
do
|
{
|
// Kahan Summation
|
// NOTE (ruegg): do NOT optimize. Now, how to tell that the compiler?
|
current = nextSummand();
|
double y = current - compensation;
|
double t = sum + y;
|
compensation = t - sum;
|
compensation -= y;
|
sum = t;
|
}
|
while (Math.Abs(sum) < Math.Abs(factor*current));
|
|
return sum;
|
}
|
|
/// <summary> Evaluates the series of Chebyshev polynomials Ti at argument x/2.
|
/// The series is given by
|
/// <pre>
|
/// N-1
|
/// - '
|
/// y = > coef[i] T (x/2)
|
/// - i
|
/// i=0
|
/// </pre>
|
/// Coefficients are stored in reverse order, i.e. the zero
|
/// order term is last in the array. Note N is the number of
|
/// coefficients, not the order.
|
/// <p/>
|
/// If coefficients are for the interval a to b, x must
|
/// have been transformed to x -> 2(2x - b - a)/(b-a) before
|
/// entering the routine. This maps x from (a, b) to (-1, 1),
|
/// over which the Chebyshev polynomials are defined.
|
/// <p/>
|
/// If the coefficients are for the inverted interval, in
|
/// which (a, b) is mapped to (1/b, 1/a), the transformation
|
/// required is x -> 2(2ab/x - b - a)/(b-a). If b is infinity,
|
/// this becomes x -> 4a/x - 1.
|
/// <p/>
|
/// SPEED:
|
/// <p/>
|
/// Taking advantage of the recurrence properties of the
|
/// Chebyshev polynomials, the routine requires one more
|
/// addition per loop than evaluating a nested polynomial of
|
/// the same degree.
|
/// </summary>
|
/// <param name="coefficients">The coefficients of the polynomial.</param>
|
/// <param name="x">Argument to the polynomial.</param>
|
/// <remarks>
|
/// Reference: https://bpm2.svn.codeplex.com/svn/Common.Numeric/Arithmetic.cs
|
/// <p/>
|
/// Marked as Deprecated in
|
/// http://people.apache.org/~isabel/mahout_site/mahout-matrix/apidocs/org/apache/mahout/jet/math/Arithmetic.html
|
/// </remarks>
|
internal static double ChebyshevA(double[] coefficients, double x)
|
{
|
// TODO: Unify, normalize, then make public
|
double b2;
|
|
int p = 0;
|
|
double b0 = coefficients[p++];
|
double b1 = 0.0;
|
int i = coefficients.Length - 1;
|
|
do
|
{
|
b2 = b1;
|
b1 = b0;
|
b0 = x*b1 - b2 + coefficients[p++];
|
}
|
while (--i > 0);
|
|
return 0.5*(b0 - b2);
|
}
|
|
/// <summary>
|
/// Summation of Chebyshev polynomials, using the Clenshaw method with Reinsch modification.
|
/// </summary>
|
/// <param name="n">The no. of terms in the sequence.</param>
|
/// <param name="coefficients">The coefficients of the Chebyshev series, length n+1.</param>
|
/// <param name="x">The value at which the series is to be evaluated.</param>
|
/// <remarks>
|
/// ORIGINAL AUTHOR:
|
/// Dr. Allan J. MacLeod; Dept. of Mathematics and Statistics, University of Paisley; High St., PAISLEY, SCOTLAND
|
/// REFERENCES:
|
/// "An error analysis of the modified Clenshaw method for evaluating Chebyshev and Fourier series"
|
/// J. Oliver, J.I.M.A., vol. 20, 1977, pp379-391
|
/// </remarks>
|
internal static double ChebyshevSum(int n, double[] coefficients, double x)
|
{
|
// TODO: Unify, normalize, then make public
|
|
// If |x| < 0.6 use the standard Clenshaw method
|
if (Math.Abs(x) < 0.6)
|
{
|
double u0 = 0.0;
|
double u1 = 0.0;
|
double u2 = 0.0;
|
double xx = x + x;
|
|
for (int i = n; i >= 0; i--)
|
{
|
u2 = u1;
|
u1 = u0;
|
u0 = xx*u1 + coefficients[i] - u2;
|
}
|
|
return (u0 - u2)/2.0;
|
}
|
|
// If ABS ( T ) > = 0.6 use the Reinsch modification
|
// T > = 0.6 code
|
if (x > 0.0)
|
{
|
double u1 = 0.0;
|
double d1 = 0.0;
|
double d2 = 0.0;
|
double xx = (x - 0.5) - 0.5;
|
xx = xx + xx;
|
|
for (int i = n; i >= 0; i--)
|
{
|
d2 = d1;
|
double u2 = u1;
|
d1 = xx*u2 + coefficients[i] + d2;
|
u1 = d1 + u2;
|
}
|
|
return (d1 + d2)/2.0;
|
}
|
else
|
{
|
// T < = -0.6 code
|
double u1 = 0.0;
|
double d1 = 0.0;
|
double d2 = 0.0;
|
double xx = (x + 0.5) + 0.5;
|
xx = xx + xx;
|
|
for (int i = n; i >= 0; i--)
|
{
|
d2 = d1;
|
double u2 = u1;
|
d1 = xx*u2 + coefficients[i] - d2;
|
u1 = d1 - u2;
|
}
|
|
return (d1 - d2)/2.0;
|
}
|
}
|
}
|
}
|
