From 733b5d5ae2c5d625211e2985ac25728ac3f54883 Mon Sep 17 00:00:00 2001
From: DongHun Kwak <dh0128.kwak@samsung.com>
Date: Mon, 21 Mar 2016 15:45:20 +0900
Subject: Imported Upstream version 1.58.0

Change-Id: If0072143aa26874812e0db6872e1efb10a3e5e94
Signed-off-by: DongHun Kwak <dh0128.kwak@samsung.com>
---
 boost/math/special_functions/gamma.hpp | 42 +++++++++++++++++++++++++++-------
 1 file changed, 34 insertions(+), 8 deletions(-)

(limited to 'boost/math/special_functions/gamma.hpp')

diff --git a/boost/math/special_functions/gamma.hpp b/boost/math/special_functions/gamma.hpp
index b6b4c574a9..1db1e4c2d3 100644
--- a/boost/math/special_functions/gamma.hpp
+++ b/boost/math/special_functions/gamma.hpp
@@ -88,10 +88,6 @@ T sinpx(T z)
    {
       z = -z;
    }
-   else
-   {
-      sign = -sign;
-   }
    T fl = floor(z);
    T dist;
    if(is_odd(fl))
@@ -1063,8 +1059,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);
@@ -1560,6 +1579,7 @@ T tgamma_ratio_imp(T x, T y, const Policy& pol)
 template <class T, class Policy>
 T gamma_p_derivative_imp(T a, T x, const Policy& pol)
 {
+   BOOST_MATH_STD_USING
    //
    // Usual error checks first:
    //
@@ -1585,8 +1605,14 @@ T gamma_p_derivative_imp(T a, T x, const Policy& pol)
       // overflow:
       return policies::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", 0, pol);
    }
-
-   f1 /= x;
+   if(f1 == 0)
+   {
+      // Underflow in calculation, use logs instead:
+      f1 = a * log(x) - x - lgamma(a, pol) - log(x);
+      f1 = exp(f1);
+   }
+   else
+      f1 /= x;
 
    return f1;
 }
-- 
cgit v1.2.3