boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-28 19:32:08 +00:00
Code Issues Projects Releases Wiki Activity

Add a convergence condition to Stirling's approximation in multiprecision tgamma() and lgamma().

Browse Source
This commit is contained in:
Christopher Kormanyos
2014-01-04 20:12:51 +01:00
parent 67fb4fd303
commit 06a2fcf887
1 changed files with 11 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

11
include/boost/math/special_functions/gamma.hpp
View File

@@ -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<T>(1) / 2) * one_over_x_pow_two_n_minus_one;
const T target_epsilon_to_break_loop = sum * pow(T(0.1F), static_cast<int>(static_cast<float>(std::numeric_limits<T>::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<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 in this divergent asymptotic expansion
// will not improve the result.
// Break from the loop.
break;
}
sum += term;
}