Browse Source

Minor special functions cleanup

pull/274/head
Christoph Ruegg 12 years ago
parent
commit
6cdafe62be
  1. 2
      CONTRIBUTORS.md
  2. 2
      src/Numerics/Numerics.csproj
  3. 5
      src/Numerics/SpecialFunctions/Beta.cs
  4. 4
      src/Numerics/SpecialFunctions/Erf.cs
  5. 38
      src/Numerics/SpecialFunctions/Evaluate.cs
  6. 56
      src/Numerics/SpecialFunctions/ExponentialIntegral.cs
  7. 10
      src/Numerics/SpecialFunctions/Factorial.cs
  8. 4
      src/Numerics/SpecialFunctions/Gamma.cs
  9. 4
      src/Numerics/SpecialFunctions/Harmonic.cs
  10. 10
      src/Numerics/SpecialFunctions/Logistic.cs
  11. 4
      src/Numerics/SpecialFunctions/ModifiedBessel.cs
  12. 38
      src/Numerics/SpecialFunctions/ModifiedStruve.cs
  13. 44
      src/Numerics/SpecialFunctions/Stability.cs
  14. 0
      src/Numerics/SpecialFunctions/TestFunctions.cs

2
CONTRIBUTORS.md

@ -33,6 +33,7 @@ Feel free to add a link to your personal site/blog and/or twitter handle.*
- Scott Stephens
- Superbest
- Anders Gustafsson
- Ashley Messer
- Gauthier Segay
- Patrick van der Velde
- Robin Neatherway
@ -49,6 +50,7 @@ Feel free to add a link to your personal site/blog and/or twitter handle.*
- Kosei ABE
- Martin Posch
- Paul Varkey
- Sunny Ahuwanya
- Till Hoffmann
- Tomas Petricek
- VicPara

2
src/Numerics/Numerics.csproj

@ -203,7 +203,7 @@
<Compile Include="SpecialFunctions\ModifiedStruve.cs" />
<Compile Include="SpecialFunctions\ModifiedBessel.cs" />
<Compile Include="SpecialFunctions\Logistic.cs" />
<Compile Include="SpecialFunctions\Test.cs" />
<Compile Include="SpecialFunctions\TestFunctions.cs" />
<Compile Include="Statistics\ArrayStatistics.cs" />
<Compile Include="Statistics\RunningStatistics.cs" />
<Compile Include="Statistics\QuantileDefinition.cs" />

5
src/Numerics/SpecialFunctions/Beta.cs

@ -33,12 +33,13 @@
// ALGLIB 2.0.1, Sergey Bochkanov
// </contribution>
using System;
using MathNet.Numerics.Properties;
// ReSharper disable CheckNamespace
namespace MathNet.Numerics
// ReSharper restore CheckNamespace
{
using System;
using Properties;
public static partial class SpecialFunctions
{

4
src/Numerics/SpecialFunctions/Erf.cs

@ -37,12 +37,12 @@
// * double ErfInvImpl(double p, double q, double s)
// </copyright>
using System;
// ReSharper disable CheckNamespace
namespace MathNet.Numerics
// ReSharper restore CheckNamespace
{
using System;
/// <summary>
/// This partial implementation of the SpecialFunctions class contains all methods related to the error function.
/// </summary>

38
src/Numerics/SpecialFunctions/Evaluate.cs

@ -32,10 +32,10 @@
// 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.
// 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
@ -44,16 +44,16 @@
// ALGLIB 2.0.1, Sergey Bochkanov
// </contribution>
// ReSharper disable CheckNamespace
namespace MathNet.Numerics
// ReSharper restore CheckNamespace
{
using System;
using System;
#if !NOSYSNUMERICS
using Complex = System.Numerics.Complex;
using System.Numerics;
#endif
// ReSharper disable CheckNamespace
namespace MathNet.Numerics
// ReSharper restore CheckNamespace
{
/// <summary>
/// Evaluation functions, useful for function approximation.
/// </summary>
@ -140,7 +140,7 @@ namespace MathNet.Numerics
compensation -= y;
sum = t;
}
while (Math.Abs(sum) < Math.Abs(factor * current));
while (Math.Abs(sum) < Math.Abs(factor*current));
return sum;
}
@ -198,11 +198,11 @@ namespace MathNet.Numerics
{
b2 = b1;
b1 = b0;
b0 = x * b1 - b2 + coefficients[p++];
b0 = x*b1 - b2 + coefficients[p++];
}
while (--i > 0);
return 0.5 * (b0 - b2);
return 0.5*(b0 - b2);
}
/// <summary>
@ -234,10 +234,10 @@ namespace MathNet.Numerics
{
u2 = u1;
u1 = u0;
u0 = xx * u1 + coefficients[i] - u2;
u0 = xx*u1 + coefficients[i] - u2;
}
return (u0 - u2) / 2.0;
return (u0 - u2)/2.0;
}
// If ABS ( T ) > = 0.6 use the Reinsch modification
@ -254,11 +254,11 @@ namespace MathNet.Numerics
{
d2 = d1;
double u2 = u1;
d1 = xx * u2 + coefficients[i] + d2;
d1 = xx*u2 + coefficients[i] + d2;
u1 = d1 + u2;
}
return (d1 + d2) / 2.0;
return (d1 + d2)/2.0;
}
else
{
@ -273,11 +273,11 @@ namespace MathNet.Numerics
{
d2 = d1;
double u2 = u1;
d1 = xx * u2 + coefficients[i] - d2;
d1 = xx*u2 + coefficients[i] - d2;
u1 = d1 - u2;
}
return (d1 - d2) / 2.0;
return (d1 - d2)/2.0;
}
}
}

56
src/Numerics/SpecialFunctions/ExponentialIntegral.cs

@ -32,12 +32,12 @@
// Ashley Messer
// </contribution>
using System;
// ReSharper disable CheckNamespace
namespace MathNet.Numerics
// ReSharper restore CheckNamespace
{
using System;
public static partial class SpecialFunctions
{
/// <summary>
@ -53,16 +53,17 @@ namespace MathNet.Numerics
/// "Advanced mathematical methods for scientists and engineers", Bender, Carl M.; Steven A. Orszag (1978). page 253
/// </para>
/// <para>
/// for x > 1 uses continued fraction approach that is often used to compute incomplete gamma.
/// for 0 < x <= 1 uses taylor series expansion
/// for x &gt; 1 uses continued fraction approach that is often used to compute incomplete gamma.
/// for 0 &lt; x &lt;= 1 uses Taylor series expansion
/// </para>
/// <para>Our unit tests suggest that the accuracy of the Exponential Integral function is correct up to 13 floating point digits.</para>
/// </remarks>
public static double ExponentialIntegral(double x, int n)
{
//parameter validation
if (n < 0 || x < 0.0 ) {
throw new ArgumentOutOfRangeException(string.Format("x and n must be positive: x={0}, n={1}", x, n));
if (n < 0 || x < 0.0)
{
throw new ArgumentOutOfRangeException(string.Format("x and n must be positive: x={0}, n={1}", x, n));
}
const double epsilon = 0.00000000000000001;
@ -79,57 +80,60 @@ namespace MathNet.Numerics
//special cases
if (n == 0)
{
result = Math.Exp( -1.0d * x ) / x;
result = Math.Exp(-1.0d*x)/x;
return result;
}
else if (x == 0.0d)
{
result = 1.0d / (ndbl - 1.0d);
result = 1.0d/(ndbl - 1.0d);
return result;
}
//general cases
//continued fraction for large x
if (x > 1.0d)
{
if (x > 1.0d)
{
b = x + ((double)n);
c = 1.0d / nearDoubleMin;
d = 1.0d / b;
c = 1.0d/nearDoubleMin;
d = 1.0d/b;
h = d;
for (i = 1; i <= maxIterations; i++)
{
a = -1.0d * ((double)i) * ((ndbl - 1.0d) + (double)i);
a = -1.0d*((double)i)*((ndbl - 1.0d) + (double)i);
b += 2.0d;
d = 1.0d / (a * d + b);
c = b + a / c;
del = c * d;
h = h * del;
d = 1.0d/(a*d + b);
c = b + a/c;
del = c*d;
h = h*del;
if (Math.Abs(del - 1.0d) < epsilon)
{
result = h * Math.Exp( -x );
result = h*Math.Exp(-x);
return result;
}
}
throw new ArithmeticException(string.Format("continued fraction failed to converge for x={0}, n={1})", x, n));
}
//series computation for small x
//series computation for small x
else
{
result = ((ndbl - 1.0d) != 0 ? 1.0 / (ndbl - 1.0d) : (-1.0d * Math.Log(x) - Constants.EulerMascheroni)); //Set first term.
result = ((ndbl - 1.0d) != 0 ? 1.0/(ndbl - 1.0d) : (-1.0d*Math.Log(x) - Constants.EulerMascheroni)); //Set first term.
for (i = 1; i <= maxIterations; i++)
{
factorial *= (-1.0d * x / ((double)i));
if (i != (ndbl - 1.0d)) { del = -factorial / (i - (ndbl - 1.0d)); }
factorial *= (-1.0d*x/((double)i));
if (i != (ndbl - 1.0d))
{
del = -factorial/(i - (ndbl - 1.0d));
}
else
{
psi = -1.0d * Constants.EulerMascheroni;
psi = -1.0d*Constants.EulerMascheroni;
for (ii = 1; ii <= (ndbl - 1.0d); ii++)
{
psi += (1.0d / ((double)ii));
psi += (1.0d/((double)ii));
}
del = factorial * (-1.0d * Math.Log(x) + psi);
del = factorial*(-1.0d*Math.Log(x) + psi);
}
result += del;
if (Math.Abs(del) < Math.Abs(result) * epsilon)
if (Math.Abs(del) < Math.Abs(result)*epsilon)
{
return result;
}

10
src/Numerics/SpecialFunctions/Factorial.cs

@ -41,8 +41,8 @@ namespace MathNet.Numerics
{
public partial class SpecialFunctions
{
private const int FactorialMaxArgument = 170;
private static double[] _factorialCache;
const int FactorialMaxArgument = 170;
static double[] _factorialCache;
/// <summary>
/// Initializes static members of the SpecialFunctions class.
@ -52,13 +52,13 @@ namespace MathNet.Numerics
InitializeFactorial();
}
private static void InitializeFactorial()
static void InitializeFactorial()
{
_factorialCache = new double[FactorialMaxArgument + 1];
_factorialCache[0] = 1.0;
for (int i = 1; i < _factorialCache.Length; i++)
{
_factorialCache[i] = _factorialCache[i - 1] * i;
_factorialCache[i] = _factorialCache[i - 1]*i;
}
}
@ -214,4 +214,4 @@ namespace MathNet.Numerics
return Math.Floor(0.5 + Math.Exp(ret));
}
}
}
}

4
src/Numerics/SpecialFunctions/Gamma.cs

@ -33,12 +33,12 @@
// ALGLIB 2.0.1, Sergey Bochkanov
// </contribution>
using System;
// ReSharper disable CheckNamespace
namespace MathNet.Numerics
// ReSharper restore CheckNamespace
{
using System;
public static partial class SpecialFunctions
{
/// <summary>

4
src/Numerics/SpecialFunctions/Harmonic.cs

@ -33,12 +33,12 @@
// ALGLIB 2.0.1, Sergey Bochkanov
// </contribution>
using System;
// ReSharper disable CheckNamespace
namespace MathNet.Numerics
// ReSharper restore CheckNamespace
{
using System;
/// <summary>
/// This partial implementation of the SpecialFunctions class contains all methods related to the harmonic function.
/// </summary>

10
src/Numerics/SpecialFunctions/Logistic.cs

@ -33,13 +33,13 @@
// ALGLIB 2.0.1, Sergey Bochkanov
// </contribution>
using System;
using MathNet.Numerics.Properties;
// ReSharper disable CheckNamespace
namespace MathNet.Numerics
// ReSharper restore CheckNamespace
{
using System;
using Properties;
/// <summary>
/// This partial implementation of the SpecialFunctions class contains all methods related to the logistic function.
/// </summary>
@ -52,7 +52,7 @@ namespace MathNet.Numerics
/// <returns>The logistic function of <paramref name="p"/>.</returns>
public static double Logistic(double p)
{
return 1.0 / (Math.Exp(-p) + 1.0);
return 1.0/(Math.Exp(-p) + 1.0);
}
/// <summary>
@ -68,7 +68,7 @@ namespace MathNet.Numerics
throw new ArgumentOutOfRangeException("p", Resources.ArgumentBetween0And1);
}
return Math.Log(p / (1.0 - p));
return Math.Log(p/(1.0 - p));
}
}
}

4
src/Numerics/SpecialFunctions/ModifiedBessel.cs

@ -44,12 +44,12 @@
// ALGLIB 2.0.1, Sergey Bochkanov
// </contribution>
using System;
// ReSharper disable CheckNamespace
namespace MathNet.Numerics
// ReSharper restore CheckNamespace
{
using System;
/// <summary>
/// This partial implementation of the SpecialFunctions class contains all methods related to the modified bessel function.
/// </summary>

38
src/Numerics/SpecialFunctions/ModifiedStruve.cs

@ -44,12 +44,12 @@
// ALGLIB 2.0.1, Sergey Bochkanov
// </contribution>
using System;
// ReSharper disable CheckNamespace
namespace MathNet.Numerics
// ReSharper restore CheckNamespace
{
using System;
/// <summary>
/// This partial implementation of the SpecialFunctions class contains all methods related to the modified bessel function.
/// </summary>
@ -246,11 +246,11 @@ namespace MathNet.Numerics
{
if (x < xlow)
{
return Constants.TwoInvPi * x;
return Constants.TwoInvPi*x;
}
double T = (4.0 * x - 24.0) / (x + 24.0);
return Constants.TwoInvPi * x * Evaluate.ChebyshevSum(nterm1, ARL0, T) * Math.Exp(x);
double T = (4.0*x - 24.0)/(x + 24.0);
return Constants.TwoInvPi*x*Evaluate.ChebyshevSum(nterm1, ARL0, T)*Math.Exp(x);
}
// Code for |xvalue| > 16
@ -261,7 +261,7 @@ namespace MathNet.Numerics
}
else
{
double T = (x - 28.0) / (4.0 - x);
double T = (x - 28.0)/(4.0 - x);
ch1 = Evaluate.ChebyshevSum(nterm2, ARL0AS, T);
}
@ -272,18 +272,18 @@ namespace MathNet.Numerics
}
else
{
double xsq = x * x;
double T = (800.0 - xsq) / (288.0 + xsq);
double xsq = x*x;
double T = (800.0 - xsq)/(288.0 + xsq);
ch2 = Evaluate.ChebyshevSum(nterm3, AI0ML0, T);
}
double test = Math.Log(ch1) - Constants.LogSqrt2Pi - Math.Log(x) / 2.0 + x;
double test = Math.Log(ch1) - Constants.LogSqrt2Pi - Math.Log(x)/2.0 + x;
if (test > Math.Log(xmax))
{
throw new ArithmeticException("ERROR IN MISCFUN FUNCTION STRVL0: ARGUMENT CAUSES OVERFLOW");
}
return Math.Exp(test) - Constants.TwoInvPi * ch2 / x;
return Math.Exp(test) - Constants.TwoInvPi*ch2/x;
}
/// <summary>
@ -488,14 +488,14 @@ namespace MathNet.Numerics
return 0.0;
}
double xsq = x * x;
double xsq = x*x;
if (x < xlow1)
{
return xsq / Constants.Pi3Over2;
return xsq/Constants.Pi3Over2;
}
double t = (4.0 * x - 24.0) / (x + 24.0);
return xsq * Evaluate.ChebyshevSum(nterm1, ARL1, t) * Math.Exp(x) / Constants.Pi3Over2;
double t = (4.0*x - 24.0)/(x + 24.0);
return xsq*Evaluate.ChebyshevSum(nterm1, ARL1, t)*Math.Exp(x)/Constants.Pi3Over2;
}
// CODE FOR |x| > 16
@ -506,7 +506,7 @@ namespace MathNet.Numerics
}
else
{
double t = (x - 30.0) / (2.0 - x);
double t = (x - 30.0)/(2.0 - x);
ch1 = Evaluate.ChebyshevSum(nterm2, ARL1AS, t);
}
@ -517,18 +517,18 @@ namespace MathNet.Numerics
}
else
{
double xsq = x * x;
double t = (800.0 - xsq) / (288.0 + xsq);
double xsq = x*x;
double t = (800.0 - xsq)/(288.0 + xsq);
ch2 = Evaluate.ChebyshevSum(nterm3, AI1ML1, t);
}
double test = Math.Log(ch1) - Constants.LogSqrt2Pi - Math.Log(x) / 2.0 + x;
double test = Math.Log(ch1) - Constants.LogSqrt2Pi - Math.Log(x)/2.0 + x;
if (test > Math.Log(xmax))
{
throw new ArithmeticException("ERROR IN MISCFUN FUNCTION STRVL1: ARGUMENT CAUSES OVERFLOW");
}
return Math.Exp(test) - Constants.TwoInvPi * ch2;
return Math.Exp(test) - Constants.TwoInvPi*ch2;
}
/// <summary>

44
src/Numerics/SpecialFunctions/Stability.cs

@ -28,16 +28,16 @@
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
// ReSharper disable CheckNamespace
namespace MathNet.Numerics
// ReSharper restore CheckNamespace
{
using System;
using System;
#if !NOSYSNUMERICS
using Complex = System.Numerics.Complex;
using System.Numerics;
#endif
// ReSharper disable CheckNamespace
namespace MathNet.Numerics
// ReSharper restore CheckNamespace
{
public partial class SpecialFunctions
{
/// <summary>
@ -81,15 +81,15 @@ namespace MathNet.Numerics
{
if (a.Magnitude > b.Magnitude)
{
var r = b.Magnitude / a.Magnitude;
return a.Magnitude * Math.Sqrt(1 + (r * r));
var r = b.Magnitude/a.Magnitude;
return a.Magnitude*Math.Sqrt(1 + (r*r));
}
if (b != 0.0)
{
// NOTE (ruegg): not "!b.AlmostZero()" to avoid convergence issues (e.g. in SVD algorithm)
var r = a.Magnitude / b.Magnitude;
return b.Magnitude * Math.Sqrt(1 + (r * r));
var r = a.Magnitude/b.Magnitude;
return b.Magnitude*Math.Sqrt(1 + (r*r));
}
return 0d;
@ -105,15 +105,15 @@ namespace MathNet.Numerics
{
if (a.Magnitude > b.Magnitude)
{
var r = b.Magnitude / a.Magnitude;
return a.Magnitude * (float)Math.Sqrt(1 + (r * r));
var r = b.Magnitude/a.Magnitude;
return a.Magnitude*(float)Math.Sqrt(1 + (r*r));
}
if (b != 0.0f)
{
// NOTE (ruegg): not "!b.AlmostZero()" to avoid convergence issues (e.g. in SVD algorithm)
var r = a.Magnitude / b.Magnitude;
return b.Magnitude * (float)Math.Sqrt(1 + (r * r));
var r = a.Magnitude/b.Magnitude;
return b.Magnitude*(float)Math.Sqrt(1 + (r*r));
}
return 0f;
@ -129,15 +129,15 @@ namespace MathNet.Numerics
{
if (Math.Abs(a) > Math.Abs(b))
{
double r = b / a;
return Math.Abs(a) * Math.Sqrt(1 + (r * r));
double r = b/a;
return Math.Abs(a)*Math.Sqrt(1 + (r*r));
}
if (b != 0.0)
{
// NOTE (ruegg): not "!b.AlmostZero()" to avoid convergence issues (e.g. in SVD algorithm)
double r = a / b;
return Math.Abs(b) * Math.Sqrt(1 + (r * r));
double r = a/b;
return Math.Abs(b)*Math.Sqrt(1 + (r*r));
}
return 0d;
@ -153,15 +153,15 @@ namespace MathNet.Numerics
{
if (Math.Abs(a) > Math.Abs(b))
{
float r = b / a;
return Math.Abs(a) * (float)Math.Sqrt(1 + (r * r));
float r = b/a;
return Math.Abs(a)*(float)Math.Sqrt(1 + (r*r));
}
if (b != 0.0)
{
// NOTE (ruegg): not "!b.AlmostZero()" to avoid convergence issues (e.g. in SVD algorithm)
float r = a / b;
return Math.Abs(b) * (float)Math.Sqrt(1 + (r * r));
float r = a/b;
return Math.Abs(b)*(float)Math.Sqrt(1 + (r*r));
}
return 0f;

0
src/Numerics/SpecialFunctions/Test.cs → src/Numerics/SpecialFunctions/TestFunctions.cs

Loading…
Cancel
Save