From 25a71ff4b4611725d32057f9c1bf423ab2c938fd Mon Sep 17 00:00:00 2001
From: Ashley Messer <glyphard@gmail.com>
Date: Sat, 3 Jan 2015 13:04:24 -0800
Subject: [PATCH] added correct handling for a=0 for regularized incomplete
 gamma function, upper and lower for issue #228 also added checksum test

---
 src/Numerics/SpecialFunctions/Gamma.cs        | 10 ++-----
 .../SpecialFunctionsTests/GammaTests.cs       | 29 +++++++++++++++++++
 2 files changed, 32 insertions(+), 7 deletions(-)

diff --git a/src/Numerics/SpecialFunctions/Gamma.cs b/src/Numerics/SpecialFunctions/Gamma.cs
index 043c8d5f..9fde572a 100644
--- a/src/Numerics/SpecialFunctions/Gamma.cs
+++ b/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;
diff --git a/src/UnitTests/SpecialFunctionsTests/GammaTests.cs b/src/UnitTests/SpecialFunctionsTests/GammaTests.cs
index 823f5cb4..fa3a3a28 100644
--- a/src/UnitTests/SpecialFunctionsTests/GammaTests.cs
+++ b/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);
+        }
+
     }
 }