From 5af7e8fdb80c203d8988a7a8d2d852c8f6d2aef7 Mon Sep 17 00:00:00 2001
From: jzmaddock <jzmaddock@gmail.com>
Date: Sat, 20 Dec 2014 19:18:25 +0000
Subject: [PATCH] [incomplete gamma] Fix corner cases identified by Rocco
 Romeo.

---
 .../boost/math/special_functions/gamma.hpp    | 27 +++++++++++++++++--
 test/test_igamma.hpp                          |  5 ++++
 2 files changed, 30 insertions(+), 2 deletions(-)

diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp
index b6b4c574a..7de524f2d 100644
--- a/include/boost/math/special_functions/gamma.hpp
+++ b/include/boost/math/special_functions/gamma.hpp
@@ -1063,8 +1063,31 @@ T gamma_incomplete_imp(T a, T x, bool normalised, bool invert,
       {
          result = gamma_incomplete_imp(a, x, true, invert, pol, p_derivative);
          if(result == 0)
-            return policies::raise_evaluation_error<T>(function, "Obtained %1% for the incomplete gamma function, but in truth we don't really know what the answer is...", result, pol);
-         result = log(result) + boost::math::lgamma(a, pol);
+         {
+            if(invert)
+            {
+               // Try http://functions.wolfram.com/06.06.06.0039.01
+               result = 1 + 1 / (12 * a) + 1 / (288 * a * a);
+               result = log(result) - a + (a - 0.5f) * log(a) + log(boost::math::constants::root_two_pi<T>());
+               if(p_derivative)
+                  *p_derivative = exp(a * log(x) - x);
+            }
+            else
+            {
+               // This is method 2 below, done in logs, we're really outside the
+               // range of this method, but since the result is almost certainly
+               // infinite, we should probably be OK:
+               result = a * log(x) - x;
+               if(p_derivative)
+                  *p_derivative = exp(result);
+               T init_value = 0;
+               result += log(detail::lower_gamma_series(a, x, pol, init_value) / a);
+            }
+         }
+         else
+         {
+            result = log(result) + boost::math::lgamma(a, pol);
+         }
       }
       if(result > tools::log_max_value<T>())
          return policies::raise_overflow_error<T>(function, 0, pol);
diff --git a/test/test_igamma.hpp b/test/test_igamma.hpp
index 40b59d0d6..65dd539a7 100644
--- a/test/test_igamma.hpp
+++ b/test/test_igamma.hpp
@@ -215,6 +215,11 @@ void test_spots(T)
    {
       BOOST_CHECK_EQUAL(::boost::math::tgamma(static_cast<T>(176), static_cast<T>(100)), std::numeric_limits<T>::infinity());
       //BOOST_CHECK_THROW(::boost::math::tgamma(static_cast<T>(176), static_cast<T>(100), throw_policy()), std::overflow_error);
+      BOOST_CHECK_EQUAL(::boost::math::tgamma(static_cast<T>(530), static_cast<T>(2000)), std::numeric_limits<T>::infinity());
+      BOOST_CHECK_EQUAL(::boost::math::tgamma(static_cast<T>(740), static_cast<T>(2500)), std::numeric_limits<T>::infinity());
+      BOOST_CHECK_EQUAL(::boost::math::tgamma(static_cast<T>(530.5), static_cast<T>(2000)), std::numeric_limits<T>::infinity());
+      BOOST_CHECK_EQUAL(::boost::math::tgamma(static_cast<T>(740.5), static_cast<T>(2500)), std::numeric_limits<T>::infinity());
+      BOOST_CHECK_EQUAL(::boost::math::tgamma_lower(static_cast<T>(10000.0f), static_cast<T>(10000.0f / 4)), std::numeric_limits<T>::infinity());
    }
    if(std::numeric_limits<T>::max_exponent >= 1024)
    {