diff options
Diffstat (limited to 'boost/math/special_functions/gamma.hpp')
| -rw-r--r-- | boost/math/special_functions/gamma.hpp | 42 |
1 files changed, 34 insertions, 8 deletions
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; } |