Browse Source

added correct handling for a=0 for regularized incomplete gamma function, upper and lower for issue #228

also added checksum test
cuda
Ashley Messer 12 years ago
parent
commit
25a71ff4b4
  1. 10
      src/Numerics/SpecialFunctions/Gamma.cs
  2. 29
      src/UnitTests/SpecialFunctionsTests/GammaTests.cs

10
src/Numerics/SpecialFunctions/Gamma.cs

@ -165,12 +165,8 @@ namespace MathNet.Numerics
const double big = 4503599627370496.0;
const double bigInv = 2.22044604925031308085e-16;
if (x <= 0d || a <= 0d)
{
return 1d;
}
if (x < 1d || x < a)
if (x < 1d || x <= a)
{
return 1d - GammaLowerRegularized(a, x);
}
@ -279,8 +275,8 @@ namespace MathNet.Numerics
{
if (x.AlmostEqual(0.0))
{
// either 0 or 1, depending on the limit direction
return double.NaN;
//use right hand limit value because so that regularized upper/lower gamma definition holds.
return 1d;
}
return 1d;

29
src/UnitTests/SpecialFunctionsTests/GammaTests.cs

@ -233,5 +233,34 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests
{
AssertHelpers.AlmostEqualRelative(f, SpecialFunctions.GammaUpperIncomplete(a, x), 13);
}
[TestCase(0.0, 0.0, 0.0)]
[TestCase(0.0, 1.0, 0.0)]
[TestCase(0.0, 2.0, 0.0)]
public void GammaUpperIncomplete_SpecialCases(double a, double x, double f)
{
double actual = SpecialFunctions.GammaUpperRegularized(a, x);
AssertHelpers.AlmostEqualRelative(f, actual, 13);
}
[TestCase(0.0, 0.0, 1.0)]
[TestCase(0.0, 1.0, 1.0)]
[TestCase(0.0, 2.0, 1.0)]
public void GammaLowerIncomplete_SpecialCases(double a, double x, double f)
{
double actual = SpecialFunctions.GammaLowerRegularized(a, x);
AssertHelpers.AlmostEqualRelative(f, actual, 13);
}
[TestCase(0.0, 0.0, 1.0)]
[TestCase(0.0, 1.0, 1.0)]
[TestCase(0.0, 2.0, 1.0)]
public void GammaIncompleteRegularized_CheckSum(double a, double x, double f)
{
double actualLower = SpecialFunctions.GammaLowerRegularized(a, x);
double actualUpper = SpecialFunctions.GammaUpperRegularized(a, x);
AssertHelpers.AlmostEqualRelative(f, actualLower + actualUpper, 13);
}
}
}

Loading…
Cancel
Save