summaryrefslogtreecommitdiff
path: root/boost/math/special_functions/gamma.hpp
diff options
context:
space:
mode:
authorChanho Park <chanho61.park@samsung.com>2014-12-11 09:55:56 (GMT)
committerChanho Park <chanho61.park@samsung.com>2014-12-11 09:55:56 (GMT)
commit08c1e93fa36a49f49325a07fe91ff92c964c2b6c (patch)
tree7a7053ceb8874b28ec4b868d4c49b500008a102e /boost/math/special_functions/gamma.hpp
parentbb4dd8289b351fae6b55e303f189127a394a1edd (diff)
downloadboost-08c1e93fa36a49f49325a07fe91ff92c964c2b6c.zip
boost-08c1e93fa36a49f49325a07fe91ff92c964c2b6c.tar.gz
boost-08c1e93fa36a49f49325a07fe91ff92c964c2b6c.tar.bz2
Imported Upstream version 1.57.0upstream/1.57.0
Diffstat (limited to 'boost/math/special_functions/gamma.hpp')
-rw-r--r--boost/math/special_functions/gamma.hpp562
1 files changed, 456 insertions, 106 deletions
diff --git a/boost/math/special_functions/gamma.hpp b/boost/math/special_functions/gamma.hpp
index 86d15b7..b6b4c57 100644
--- a/boost/math/special_functions/gamma.hpp
+++ b/boost/math/special_functions/gamma.hpp
@@ -1,6 +1,8 @@
-// Copyright John Maddock 2006-7.
-// Copyright Paul A. Bristow 2007.
+// Copyright John Maddock 2006-7, 2013-14.
+// Copyright Paul A. Bristow 2007, 2013-14.
+// Copyright Nikhar Agrawal 2013-14
+// Copyright Christopher Kormanyos 2013-14
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0. (See accompanying file
@@ -14,17 +16,6 @@
#endif
#include <boost/config.hpp>
-#ifdef BOOST_MSVC
-# pragma warning(push)
-# pragma warning(disable: 4127 4701)
-// // For lexical_cast, until fixed in 1.35?
-// // conditional expression is constant &
-// // Potentially uninitialized local variable 'name' used
-#endif
-#include <boost/lexical_cast.hpp>
-#ifdef BOOST_MSVC
-# pragma warning(pop)
-#endif
#include <boost/math/tools/series.hpp>
#include <boost/math/tools/fraction.hpp>
#include <boost/math/tools/precision.hpp>
@@ -41,6 +32,7 @@
#include <boost/math/special_functions/detail/igamma_large.hpp>
#include <boost/math/special_functions/detail/unchecked_factorial.hpp>
#include <boost/math/special_functions/detail/lgamma_small.hpp>
+#include <boost/math/special_functions/bernoulli.hpp>
#include <boost/type_traits/is_convertible.hpp>
#include <boost/assert.hpp>
#include <boost/mpl/greater.hpp>
@@ -50,12 +42,6 @@
#include <boost/config/no_tr1/cmath.hpp>
#include <algorithm>
-#ifdef BOOST_MATH_INSTRUMENT
-#include <iostream>
-#include <iomanip>
-#include <typeinfo>
-#endif
-
#ifdef BOOST_MSVC
# pragma warning(push)
# pragma warning(disable: 4702) // unreachable code (return after domain_error throw).
@@ -153,7 +139,7 @@ T gamma_imp(T z, const Policy& pol, const Lanczos& l)
result = gamma_imp(T(-z), pol, l) * sinpx(z);
BOOST_MATH_INSTRUMENT_VARIABLE(result);
if((fabs(result) < 1) && (tools::max_value<T>() * fabs(result) < boost::math::constants::pi<T>()))
- return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
+ return -boost::math::sign(result) * policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
result = -boost::math::constants::pi<T>() / result;
if(result == 0)
return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);
@@ -176,30 +162,36 @@ T gamma_imp(T z, const Policy& pol, const Lanczos& l)
result *= unchecked_factorial<T>(itrunc(z, pol) - 1);
BOOST_MATH_INSTRUMENT_VARIABLE(result);
}
+ else if (z < tools::root_epsilon<T>())
+ {
+ if (z < 1 / tools::max_value<T>())
+ result = policies::raise_overflow_error<T>(function, 0, pol);
+ result *= 1 / z - constants::euler<T>();
+ }
else
{
result *= Lanczos::lanczos_sum(z);
+ T zgh = (z + static_cast<T>(Lanczos::g()) - boost::math::constants::half<T>());
+ T lzgh = log(zgh);
BOOST_MATH_INSTRUMENT_VARIABLE(result);
BOOST_MATH_INSTRUMENT_VARIABLE(tools::log_max_value<T>());
- if(z * log(z) > tools::log_max_value<T>())
+ if(z * lzgh > tools::log_max_value<T>())
{
// we're going to overflow unless this is done with care:
- T zgh = (z + static_cast<T>(Lanczos::g()) - boost::math::constants::half<T>());
BOOST_MATH_INSTRUMENT_VARIABLE(zgh);
- if(log(zgh) * z / 2 > tools::log_max_value<T>())
- return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
+ if(lzgh * z / 2 > tools::log_max_value<T>())
+ return boost::math::sign(result) * policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
T hp = pow(zgh, (z / 2) - T(0.25));
BOOST_MATH_INSTRUMENT_VARIABLE(hp);
result *= hp / exp(zgh);
BOOST_MATH_INSTRUMENT_VARIABLE(result);
if(tools::max_value<T>() / hp < result)
- return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
+ return boost::math::sign(result) * policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
result *= hp;
BOOST_MATH_INSTRUMENT_VARIABLE(result);
}
else
{
- T zgh = (z + static_cast<T>(Lanczos::g()) - boost::math::constants::half<T>());
BOOST_MATH_INSTRUMENT_VARIABLE(zgh);
BOOST_MATH_INSTRUMENT_VARIABLE(pow(zgh, z - boost::math::constants::half<T>()));
BOOST_MATH_INSTRUMENT_VARIABLE(exp(zgh));
@@ -230,7 +222,7 @@ T lgamma_imp(T z, const Policy& pol, const Lanczos& l, int* sign = 0)
T result = 0;
int sresult = 1;
- if(z <= 0)
+ if(z <= -tools::root_epsilon<T>())
{
// reflection formula:
if(floor(z) == z)
@@ -248,6 +240,17 @@ T lgamma_imp(T z, const Policy& pol, const Lanczos& l, int* sign = 0)
}
result = log(boost::math::constants::pi<T>()) - lgamma_imp(z, pol, l) - log(t);
}
+ else if (z < tools::root_epsilon<T>())
+ {
+ if (0 == z)
+ return policies::raise_pole_error<T>(function, "Evaluation of lgamma at %1%.", z, pol);
+ if (fabs(z) < 1 / tools::max_value<T>())
+ result = -log(fabs(z));
+ else
+ result = log(fabs(1 / z - constants::euler<T>()));
+ if (z < 0)
+ sresult = -1;
+ }
else if(z < 15)
{
typedef typename policies::precision<T, Policy>::type precision_type;
@@ -266,7 +269,7 @@ T lgamma_imp(T z, const Policy& pol, const Lanczos& l, int* sign = 0)
>::type tag_type;
result = lgamma_small_imp<T>(z, T(z - 1), T(z - 2), tag_type(), pol, l);
}
- else if((z >= 3) && (z < 100))
+ else if((z >= 3) && (z < 100) && (std::numeric_limits<T>::max_exponent >= 1024))
{
// taking the log of tgamma reduces the error, no danger of overflow here:
result = log(gamma_imp(z, pol, l));
@@ -353,96 +356,271 @@ inline T lower_gamma_series(T a, T z, const Policy& pol, T init_value = 0)
}
//
-// Fully generic tgamma and lgamma use the incomplete partial
-// sums added together:
+// Fully generic tgamma and lgamma use Stirling's approximation
+// with Bernoulli numbers.
//
+template<class T>
+std::size_t highest_bernoulli_index()
+{
+ const float digits10_of_type = (std::numeric_limits<T>::is_specialized
+ ? static_cast<float>(std::numeric_limits<T>::digits10)
+ : static_cast<float>(boost::math::tools::digits<T>() * 0.301F));
+
+ // Find the high index n for Bn to produce the desired precision in Stirling's calculation.
+ return static_cast<std::size_t>(18.0F + (0.6F * digits10_of_type));
+}
+
+template<class T>
+T minimum_argument_for_bernoulli_recursion()
+{
+ const float digits10_of_type = (std::numeric_limits<T>::is_specialized
+ ? static_cast<float>(std::numeric_limits<T>::digits10)
+ : static_cast<float>(boost::math::tools::digits<T>() * 0.301F));
+
+ return T(digits10_of_type * 1.7F);
+}
+
+// Forward declaration of the lgamma_imp template specialization.
template <class T, class Policy>
-T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos& l)
+T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sign = 0);
+
+template <class T, class Policy>
+T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&)
{
- static const char* function = "boost::math::tgamma<%1%>(%1%)";
BOOST_MATH_STD_USING
- if((z <= 0) && (floor(z) == z))
- return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
- if(z <= -20)
+
+ static const char* function = "boost::math::tgamma<%1%>(%1%)";
+
+ // Check if the argument of tgamma is identically zero.
+ const bool is_at_zero = (z == 0);
+
+ if(is_at_zero)
+ return policies::raise_domain_error<T>(function, "Evaluation of tgamma at zero %1%.", z, pol);
+
+ const bool b_neg = (z < 0);
+
+ const bool floor_of_z_is_equal_to_z = (floor(z) == z);
+
+ // Special case handling of small factorials:
+ if((!b_neg) && floor_of_z_is_equal_to_z && (z < boost::math::max_factorial<T>::value))
{
- T result = gamma_imp(T(-z), pol, l) * sinpx(z);
- if((fabs(result) < 1) && (tools::max_value<T>() * fabs(result) < boost::math::constants::pi<T>()))
- return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
- result = -boost::math::constants::pi<T>() / result;
- if(result == 0)
- return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);
- if((boost::math::fpclassify)(result) == (int)FP_SUBNORMAL)
- return policies::raise_denorm_error<T>(function, "Result of tgamma is denormalized.", result, pol);
- return result;
+ return boost::math::unchecked_factorial<T>(itrunc(z) - 1);
}
- //
- // The upper gamma fraction is *very* slow for z < 6, actually it's very
- // slow to converge everywhere but recursing until z > 6 gets rid of the
- // worst of it's behaviour.
- //
- T prefix = 1;
- while(z < 6)
+
+ // Make a local, unsigned copy of the input argument.
+ T zz((!b_neg) ? z : -z);
+
+ // Special case for ultra-small z:
+ if(zz < tools::cbrt_epsilon<T>())
{
- prefix /= z;
- z += 1;
+ const T a0(1);
+ const T a1(boost::math::constants::euler<T>());
+ const T six_euler_squared((boost::math::constants::euler<T>() * boost::math::constants::euler<T>()) * 6);
+ const T a2((six_euler_squared - boost::math::constants::pi_sqr<T>()) / 12);
+
+ const T inverse_tgamma_series = z * ((a2 * z + a1) * z + a0);
+
+ return 1 / inverse_tgamma_series;
}
- BOOST_MATH_INSTRUMENT_CODE(prefix);
- if((floor(z) == z) && (z < max_factorial<T>::value))
+
+ // Scale the argument up for the calculation of lgamma,
+ // and use downward recursion later for the final result.
+ const T min_arg_for_recursion = minimum_argument_for_bernoulli_recursion<T>();
+
+ int n_recur;
+
+ if(zz < min_arg_for_recursion)
{
- prefix *= unchecked_factorial<T>(itrunc(z, pol) - 1);
+ n_recur = boost::math::itrunc(min_arg_for_recursion - zz) + 1;
+
+ zz += n_recur;
}
else
{
- prefix = prefix * pow(z / boost::math::constants::e<T>(), z);
- BOOST_MATH_INSTRUMENT_CODE(prefix);
- T sum = detail::lower_gamma_series(z, z, pol) / z;
- BOOST_MATH_INSTRUMENT_CODE(sum);
- sum += detail::upper_gamma_fraction(z, z, ::boost::math::policies::get_epsilon<T, Policy>());
- BOOST_MATH_INSTRUMENT_CODE(sum);
- if(fabs(tools::max_value<T>() / prefix) < fabs(sum))
+ n_recur = 0;
+ }
+
+ const T log_gamma_value = lgamma_imp(zz, pol, lanczos::undefined_lanczos());
+
+ if(log_gamma_value > tools::log_max_value<T>())
+ return policies::raise_overflow_error<T>(function, 0, pol);
+
+ T gamma_value = exp(log_gamma_value);
+
+ // Rescale the result using downward recursion if necessary.
+ if(n_recur)
+ {
+ // The order of divides is important, if we keep subtracting 1 from zz
+ // we DO NOT get back to z (cancellation error). Further if z < epsilon
+ // we would end up dividing by zero. Also in order to prevent spurious
+ // overflow with the first division, we must save dividing by |z| till last,
+ // so the optimal order of divides is z+1, z+2, z+3...z+n_recur-1,z.
+ zz = fabs(z) + 1;
+ for(int k = 1; k < n_recur; ++k)
+ {
+ gamma_value /= zz;
+ zz += 1;
+ }
+ gamma_value /= fabs(z);
+ }
+
+ // Return the result, accounting for possible negative arguments.
+ if(b_neg)
+ {
+ // Provide special error analysis for:
+ // * arguments in the neighborhood of a negative integer
+ // * arguments exactly equal to a negative integer.
+
+ // Check if the argument of tgamma is exactly equal to a negative integer.
+ if(floor_of_z_is_equal_to_z)
+ return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
+
+ gamma_value *= sinpx(z);
+
+ BOOST_MATH_INSTRUMENT_VARIABLE(gamma_value);
+
+ const bool result_is_too_large_to_represent = ( (abs(gamma_value) < 1)
+ && ((tools::max_value<T>() * abs(gamma_value)) < boost::math::constants::pi<T>()));
+
+ if(result_is_too_large_to_represent)
return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
- BOOST_MATH_INSTRUMENT_CODE((sum * prefix));
- return sum * prefix;
+
+ gamma_value = -boost::math::constants::pi<T>() / gamma_value;
+ BOOST_MATH_INSTRUMENT_VARIABLE(gamma_value);
+
+ if(gamma_value == 0)
+ return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);
+
+ if((boost::math::fpclassify)(gamma_value) == static_cast<int>(FP_SUBNORMAL))
+ return policies::raise_denorm_error<T>(function, "Result of tgamma is denormalized.", gamma_value, pol);
}
- return prefix;
+
+ return gamma_value;
}
template <class T, class Policy>
-T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos& l, int*sign)
+T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sign)
{
BOOST_MATH_STD_USING
static const char* function = "boost::math::lgamma<%1%>(%1%)";
- T result = 0;
- int sresult = 1;
- if(z <= 0)
+
+ // Check if the argument of lgamma is identically zero.
+ const bool is_at_zero = (z == 0);
+
+ if(is_at_zero)
+ return policies::raise_domain_error<T>(function, "Evaluation of lgamma at zero %1%.", z, pol);
+
+ const bool b_neg = (z < 0);
+
+ const bool floor_of_z_is_equal_to_z = (floor(z) == z);
+
+ // Special case handling of small factorials:
+ if((!b_neg) && floor_of_z_is_equal_to_z && (z < boost::math::max_factorial<T>::value))
{
- if(floor(z) == z)
- return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
- T t = detail::sinpx(z);
- z = -z;
+ return log(boost::math::unchecked_factorial<T>(itrunc(z) - 1));
+ }
+
+ // Make a local, unsigned copy of the input argument.
+ T zz((!b_neg) ? z : -z);
+
+ const T min_arg_for_recursion = minimum_argument_for_bernoulli_recursion<T>();
+
+ T log_gamma_value;
+
+ if (zz < min_arg_for_recursion)
+ {
+ // Here we simply take the logarithm of tgamma(). This is somewhat
+ // inefficient, but simple. The rationale is that the argument here
+ // is relatively small and overflow is not expected to be likely.
+ if (z > -tools::root_epsilon<T>())
+ {
+ // Reflection formula may fail if z is very close to zero, let the series
+ // expansion for tgamma close to zero do the work:
+ log_gamma_value = log(abs(gamma_imp(z, pol, lanczos::undefined_lanczos())));
+ if (sign)
+ {
+ *sign = z < 0 ? -1 : 1;
+ }
+ return log_gamma_value;
+ }
+ else
+ {
+ // No issue with spurious overflow in reflection formula,
+ // just fall through to regular code:
+ log_gamma_value = log(abs(gamma_imp(zz, pol, lanczos::undefined_lanczos())));
+ }
+ }
+ else
+ {
+ // Perform the Bernoulli series expansion of Stirling's approximation.
+
+ const std::size_t number_of_bernoullis_b2n = highest_bernoulli_index<T>();
+
+ T one_over_x_pow_two_n_minus_one = 1 / zz;
+ const T one_over_x2 = one_over_x_pow_two_n_minus_one * one_over_x_pow_two_n_minus_one;
+ T sum = (boost::math::bernoulli_b2n<T>(1) / 2) * one_over_x_pow_two_n_minus_one;
+ const T target_epsilon_to_break_loop = (sum * boost::math::tools::epsilon<T>()) * T(1.0E-10F);
+
+ for(std::size_t n = 2U; n < number_of_bernoullis_b2n; ++n)
+ {
+ one_over_x_pow_two_n_minus_one *= one_over_x2;
+
+ const std::size_t n2 = static_cast<std::size_t>(n * 2U);
+
+ const T term = (boost::math::bernoulli_b2n<T>(static_cast<int>(n)) * one_over_x_pow_two_n_minus_one) / (n2 * (n2 - 1U));
+
+ if((n >= 8U) && (abs(term) < target_epsilon_to_break_loop))
+ {
+ // We have reached the desired precision in Stirling's expansion.
+ // Adding additional terms to the sum of this divergent asymptotic
+ // expansion will not improve the result.
+
+ // Break from the loop.
+ break;
+ }
+
+ sum += term;
+ }
+
+ // Complete Stirling's approximation.
+ const T half_ln_two_pi = log(boost::math::constants::two_pi<T>()) / 2;
+
+ log_gamma_value = ((((zz - boost::math::constants::half<T>()) * log(zz)) - zz) + half_ln_two_pi) + sum;
+ }
+
+ int sign_of_result = 1;
+
+ if(b_neg)
+ {
+ // Provide special error analysis if the argument is exactly
+ // equal to a negative integer.
+
+ // Check if the argument of lgamma is exactly equal to a negative integer.
+ if(floor_of_z_is_equal_to_z)
+ return policies::raise_pole_error<T>(function, "Evaluation of lgamma at a negative integer %1%.", z, pol);
+
+ T t = sinpx(z);
+
if(t < 0)
{
t = -t;
}
else
{
- sresult = -sresult;
+ sign_of_result = -sign_of_result;
}
- result = log(boost::math::constants::pi<T>()) - lgamma_imp(z, pol, l, 0) - log(t);
- }
- else if((z != 1) && (z != 2))
- {
- T limit = (std::max)(T(z+1), T(10));
- T prefix = z * log(limit) - limit;
- T sum = detail::lower_gamma_series(z, limit, pol) / z;
- sum += detail::upper_gamma_fraction(z, limit, ::boost::math::policies::get_epsilon<T, Policy>());
- result = log(sum) + prefix;
+
+ log_gamma_value = - log_gamma_value
+ + log(boost::math::constants::pi<T>())
+ - log(t);
}
- if(sign)
- *sign = sresult;
- return result;
+
+ if(sign != static_cast<int*>(0U)) { *sign = sign_of_result; }
+
+ return log_gamma_value;
}
+
//
// This helper calculates tgamma(dz+1)-1 without cancellation errors,
// used by the upper incomplete gamma with z < 1:
@@ -604,7 +782,7 @@ T full_igamma_prefix(T a, T z, const Policy& pol)
// rather than before it...
//
if((boost::math::fpclassify)(prefix) == (int)FP_INFINITE)
- policies::raise_overflow_error<T>("boost::math::detail::full_igamma_prefix<%1%>(%1%, %1%)", "Result of incomplete gamma function is too large to represent.", pol);
+ return policies::raise_overflow_error<T>("boost::math::detail::full_igamma_prefix<%1%>(%1%, %1%)", "Result of incomplete gamma function is too large to represent.", pol);
return prefix;
}
@@ -842,9 +1020,9 @@ T gamma_incomplete_imp(T a, T x, bool normalised, bool invert,
{
static const char* function = "boost::math::gamma_p<%1%>(%1%, %1%)";
if(a <= 0)
- policies::raise_domain_error<T>(function, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol);
+ return policies::raise_domain_error<T>(function, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol);
if(x < 0)
- policies::raise_domain_error<T>(function, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol);
+ return policies::raise_domain_error<T>(function, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol);
BOOST_MATH_STD_USING
@@ -852,10 +1030,51 @@ T gamma_incomplete_imp(T a, T x, bool normalised, bool invert,
T result = 0; // Just to avoid warning C4701: potentially uninitialized local variable 'result' used
+ if(a >= max_factorial<T>::value && !normalised)
+ {
+ //
+ // When we're computing the non-normalized incomplete gamma
+ // and a is large the result is rather hard to compute unless
+ // we use logs. There are really two options - if x is a long
+ // way from a in value then we can reliably use methods 2 and 4
+ // below in logarithmic form and go straight to the result.
+ // Otherwise we let the regularized gamma take the strain
+ // (the result is unlikely to unerflow in the central region anyway)
+ // and combine with lgamma in the hopes that we get a finite result.
+ //
+ if(invert && (a * 4 < x))
+ {
+ // This is method 4 below, done in logs:
+ result = a * log(x) - x;
+ if(p_derivative)
+ *p_derivative = exp(result);
+ result += log(upper_gamma_fraction(a, x, policies::get_epsilon<T, Policy>()));
+ }
+ else if(!invert && (a > 4 * x))
+ {
+ // This is method 2 below, done in logs:
+ 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 = 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(result > tools::log_max_value<T>())
+ return policies::raise_overflow_error<T>(function, 0, pol);
+ return exp(result);
+ }
+
BOOST_ASSERT((p_derivative == 0) || (normalised == true));
bool is_int, is_half_int;
- bool is_small_a = (a < 30) && (a <= x + 1);
+ bool is_small_a = (a < 30) && (a <= x + 1) && (x < tools::log_max_value<T>());
if(is_small_a)
{
T fa = floor(a);
@@ -881,6 +1100,10 @@ T gamma_incomplete_imp(T a, T x, bool normalised, bool invert,
invert = !invert;
eval_method = 1;
}
+ else if((x < tools::root_epsilon<T>()) && (a > 1))
+ {
+ eval_method = 6;
+ }
else if(x < 0.5)
{
//
@@ -994,13 +1217,39 @@ T gamma_incomplete_imp(T a, T x, bool normalised, bool invert,
*p_derivative = result;
if(result != 0)
{
+ //
+ // If we're going to be inverting the result then we can
+ // reduce the number of series evaluations by quite
+ // a few iterations if we set an initial value for the
+ // series sum based on what we'll end up subtracting it from
+ // at the end.
+ // Have to be careful though that this optimization doesn't
+ // lead to spurious numberic overflow. Note that the
+ // scary/expensive overflow checks below are more often
+ // than not bypassed in practice for "sensible" input
+ // values:
+ //
T init_value = 0;
+ bool optimised_invert = false;
if(invert)
{
- init_value = -a * (normalised ? 1 : boost::math::tgamma(a, pol)) / result;
+ init_value = (normalised ? 1 : boost::math::tgamma(a, pol));
+ if(normalised || (result >= 1) || (tools::max_value<T>() * result > init_value))
+ {
+ init_value /= result;
+ if(normalised || (a < 1) || (tools::max_value<T>() / a > init_value))
+ {
+ init_value *= -a;
+ optimised_invert = true;
+ }
+ else
+ init_value = 0;
+ }
+ else
+ init_value = 0;
}
result *= detail::lower_gamma_series(a, x, pol, init_value) / a;
- if(invert)
+ if(optimised_invert)
{
invert = false;
result = -result;
@@ -1063,6 +1312,13 @@ T gamma_incomplete_imp(T a, T x, bool normalised, bool invert,
*p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
break;
}
+ case 6:
+ {
+ // x is so small that P is necessarily very small too,
+ // use http://functions.wolfram.com/GammaBetaErf/GammaRegularized/06/01/05/01/01/
+ result = !normalised ? pow(x, a) / (a) : pow(x, a) / boost::math::tgamma(a + 1, pol);
+ result *= 1 - a * x / (a + 1);
+ }
}
if(normalised && (result > 1))
@@ -1093,9 +1349,32 @@ T gamma_incomplete_imp(T a, T x, bool normalised, bool invert,
// Ratios of two gamma functions:
//
template <class T, class Policy, class Lanczos>
-T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const Lanczos&)
+T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const Lanczos& l)
{
BOOST_MATH_STD_USING
+ if(z < tools::epsilon<T>())
+ {
+ //
+ // We get spurious numeric overflow unless we're very careful, this
+ // can occur either inside Lanczos::lanczos_sum(z) or in the
+ // final combination of terms, to avoid this, split the product up
+ // into 2 (or 3) parts:
+ //
+ // G(z) / G(L) = 1 / (z * G(L)) ; z < eps, L = z + delta = delta
+ // z * G(L) = z * G(lim) * (G(L)/G(lim)) ; lim = largest factorial
+ //
+ if(boost::math::max_factorial<T>::value < delta)
+ {
+ T ratio = tgamma_delta_ratio_imp_lanczos(delta, T(boost::math::max_factorial<T>::value - delta), pol, l);
+ ratio *= z;
+ ratio *= boost::math::unchecked_factorial<T>(boost::math::max_factorial<T>::value - 1);
+ return 1 / ratio;
+ }
+ else
+ {
+ return 1 / (z * boost::math::tgamma(z + delta, pol));
+ }
+ }
T zgh = z + Lanczos::g() - constants::half<T>();
T result;
if(fabs(delta) < 10)
@@ -1106,8 +1385,9 @@ T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const Lanczos&
{
result = pow(zgh / (zgh + delta), z - constants::half<T>());
}
- result *= pow(constants::e<T>() / (zgh + delta), delta);
+ // Split the calculation up to avoid spurious overflow:
result *= Lanczos::lanczos_sum(z) / Lanczos::lanczos_sum(T(z + delta));
+ result *= pow(constants::e<T>() / (zgh + delta), delta);
return result;
}
//
@@ -1155,10 +1435,11 @@ T tgamma_delta_ratio_imp(T z, T delta, const Policy& pol)
{
BOOST_MATH_STD_USING
- if(z <= 0)
- policies::raise_domain_error<T>("boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)", "Gamma function ratios only implemented for positive arguments (got a=%1%).", z, pol);
- if(z+delta <= 0)
- policies::raise_domain_error<T>("boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)", "Gamma function ratios only implemented for positive arguments (got b=%1%).", z+delta, pol);
+ if((z <= 0) || (z + delta <= 0))
+ {
+ // This isn't very sofisticated, or accurate, but it does work:
+ return boost::math::tgamma(z, pol) / boost::math::tgamma(z + delta, pol);
+ }
if(floor(delta) == delta)
{
@@ -1208,15 +1489,84 @@ T tgamma_delta_ratio_imp(T z, T delta, const Policy& pol)
}
template <class T, class Policy>
+T tgamma_ratio_imp(T x, T y, const Policy& pol)
+{
+ BOOST_MATH_STD_USING
+
+ if((x <= 0) || (boost::math::isinf)(x))
+ return policies::raise_domain_error<T>("boost::math::tgamma_ratio<%1%>(%1%, %1%)", "Gamma function ratios only implemented for positive arguments (got a=%1%).", x, pol);
+ if((y <= 0) || (boost::math::isinf)(y))
+ return policies::raise_domain_error<T>("boost::math::tgamma_ratio<%1%>(%1%, %1%)", "Gamma function ratios only implemented for positive arguments (got b=%1%).", y, pol);
+
+ if(x <= tools::min_value<T>())
+ {
+ // Special case for denorms...Ugh.
+ T shift = ldexp(T(1), tools::digits<T>());
+ return shift * tgamma_ratio_imp(T(x * shift), y, pol);
+ }
+
+ if((x < max_factorial<T>::value) && (y < max_factorial<T>::value))
+ {
+ // Rather than subtracting values, lets just call the gamma functions directly:
+ return boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
+ }
+ T prefix = 1;
+ if(x < 1)
+ {
+ if(y < 2 * max_factorial<T>::value)
+ {
+ // We need to sidestep on x as well, otherwise we'll underflow
+ // before we get to factor in the prefix term:
+ prefix /= x;
+ x += 1;
+ while(y >= max_factorial<T>::value)
+ {
+ y -= 1;
+ prefix /= y;
+ }
+ return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
+ }
+ //
+ // result is almost certainly going to underflow to zero, try logs just in case:
+ //
+ return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol));
+ }
+ if(y < 1)
+ {
+ if(x < 2 * max_factorial<T>::value)
+ {
+ // We need to sidestep on y as well, otherwise we'll overflow
+ // before we get to factor in the prefix term:
+ prefix *= y;
+ y += 1;
+ while(x >= max_factorial<T>::value)
+ {
+ x -= 1;
+ prefix *= x;
+ }
+ return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
+ }
+ //
+ // Result will almost certainly overflow, try logs just in case:
+ //
+ return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol));
+ }
+ //
+ // Regular case, x and y both large and similar in magnitude:
+ //
+ return boost::math::tgamma_delta_ratio(x, y - x, pol);
+}
+
+template <class T, class Policy>
T gamma_p_derivative_imp(T a, T x, const Policy& pol)
{
//
// Usual error checks first:
//
if(a <= 0)
- policies::raise_domain_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol);
+ return policies::raise_domain_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol);
if(x < 0)
- policies::raise_domain_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol);
+ return policies::raise_domain_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol);
//
// Now special cases:
//
@@ -1360,7 +1710,7 @@ inline typename tools::promote_args<T1, T2>::type
BOOST_FPU_EXCEPTION_GUARD
typedef typename tools::promote_args<T1, T2>::type result_type;
typedef typename policies::evaluation<result_type, Policy>::type value_type;
- typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
+ // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
@@ -1489,7 +1839,7 @@ inline typename tools::promote_args<T1, T2>::type
BOOST_FPU_EXCEPTION_GUARD
typedef typename tools::promote_args<T1, T2>::type result_type;
typedef typename policies::evaluation<result_type, Policy>::type value_type;
- typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
+ // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
@@ -1520,7 +1870,7 @@ inline typename tools::promote_args<T1, T2>::type
BOOST_FPU_EXCEPTION_GUARD
typedef typename tools::promote_args<T1, T2>::type result_type;
typedef typename policies::evaluation<result_type, Policy>::type value_type;
- typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
+ // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
@@ -1551,7 +1901,7 @@ inline typename tools::promote_args<T1, T2>::type
BOOST_FPU_EXCEPTION_GUARD
typedef typename tools::promote_args<T1, T2>::type result_type;
typedef typename policies::evaluation<result_type, Policy>::type value_type;
- typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
+ // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
@@ -1609,7 +1959,7 @@ inline typename tools::promote_args<T1, T2>::type
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
- return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::tgamma_delta_ratio_imp(static_cast<value_type>(a), static_cast<value_type>(static_cast<value_type>(b) - static_cast<value_type>(a)), forwarding_policy()), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)");
+ return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::tgamma_ratio_imp(static_cast<value_type>(a), static_cast<value_type>(b), forwarding_policy()), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)");
}
template <class T1, class T2>
inline typename tools::promote_args<T1, T2>::type