summaryrefslogtreecommitdiff
path: root/boost/math/special_functions/gamma.hpp
diff options
context:
space:
mode:
Diffstat (limited to 'boost/math/special_functions/gamma.hpp')
-rw-r--r--boost/math/special_functions/gamma.hpp42
1 files changed, 34 insertions, 8 deletions
diff --git a/boost/math/special_functions/gamma.hpp b/boost/math/special_functions/gamma.hpp
index b6b4c57..1db1e4c 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;
}