From 06a2fcf887e2f4915d74f4fa60886eecbb3c5906 Mon Sep 17 00:00:00 2001 From: Christopher Kormanyos Date: Sat, 4 Jan 2014 20:12:51 +0100 Subject: [PATCH] Add a convergence condition to Stirling's approximation in multiprecision tgamma() and lgamma(). --- include/boost/math/special_functions/gamma.hpp | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp index dd1821df0..06bbc4641 100644 --- a/include/boost/math/special_functions/gamma.hpp +++ b/include/boost/math/special_functions/gamma.hpp @@ -480,6 +480,7 @@ T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sig 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(1) / 2) * one_over_x_pow_two_n_minus_one; + const T target_epsilon_to_break_loop = sum * pow(T(0.1F), static_cast(static_cast(std::numeric_limits::digits10) * 1.15F)); for(std::size_t n = 2U; n < number_of_bernoullis_b2n; ++n) { @@ -489,6 +490,16 @@ T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sig const T term = (boost::math::bernoulli_b2n(static_cast(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 in this divergent asymptotic expansion + // will not improve the result. + + // Break from the loop. + break; + } + sum += term; }