X-Git-Url: https://git.proxmox.com/?a=blobdiff_plain;f=ceph%2Fsrc%2Fboost%2Fboost%2Fmath%2Fspecial_functions%2Fgamma.hpp;h=bdbff31554a7ace4f047287efa0be18e0b72b79b;hb=92f5a8d42d07f9929ae4fa7e01342fe8d96808a8;hp=1903b997cbf7321c31c88aa017c2779ba1ef0d7b;hpb=a0324939f9d0e1905d5df8f57442f09dc70af83d;p=ceph.git diff --git a/ceph/src/boost/boost/math/special_functions/gamma.hpp b/ceph/src/boost/boost/math/special_functions/gamma.hpp index 1903b997c..bdbff3155 100644 --- a/ceph/src/boost/boost/math/special_functions/gamma.hpp +++ b/ceph/src/boost/boost/math/special_functions/gamma.hpp @@ -33,7 +33,7 @@ #include #include #include -#include +#include #include #include #include @@ -277,7 +277,11 @@ T lgamma_imp(T z, const Policy& pol, const Lanczos& l, int* sign = 0) T zgh = static_cast(z + Lanczos::g() - boost::math::constants::half()); result = log(zgh) - 1; result *= z - 0.5f; - result += log(Lanczos::lanczos_sum_expG_scaled(z)); + // + // Only add on the lanczos sum part if we're going to need it: + // + if(result * tools::epsilon() < 20) + result += log(Lanczos::lanczos_sum_expG_scaled(z)); } if(sign) @@ -368,13 +372,71 @@ std::size_t highest_bernoulli_index() } template -T minimum_argument_for_bernoulli_recursion() +int minimum_argument_for_bernoulli_recursion() { const float digits10_of_type = (std::numeric_limits::is_specialized ? static_cast(std::numeric_limits::digits10) : static_cast(boost::math::tools::digits() * 0.301F)); - return T(digits10_of_type * 1.7F); + const float limit = std::ceil(std::pow(1.0f / std::ldexp(1.0f, 1-boost::math::tools::digits()), 1.0f / 20.0f)); + + return (int)((std::min)(digits10_of_type * 1.7F, limit)); +} + +template +T scaled_tgamma_no_lanczos(const T& z, const Policy& pol, bool islog = false) +{ + BOOST_MATH_STD_USING + // + // Calculates tgamma(z) / (z/e)^z + // Requires that our argument is large enough for Sterling's approximation to hold. + // Used internally when combining gamma's of similar magnitude without logarithms. + // + BOOST_ASSERT(minimum_argument_for_bernoulli_recursion() <= z); + + // Perform the Bernoulli series expansion of Stirling's approximation. + + const std::size_t number_of_bernoullis_b2n = policies::get_max_series_iterations(); + + T one_over_x_pow_two_n_minus_one = 1 / z; + 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 * boost::math::tools::epsilon(); + const T half_ln_two_pi_over_z = sqrt(boost::math::constants::two_pi() / z); + T last_term = 2 * sum; + + for (std::size_t n = 2U;; ++n) + { + one_over_x_pow_two_n_minus_one *= one_over_x2; + + const std::size_t n2 = static_cast(n * 2U); + + const T term = (boost::math::bernoulli_b2n(static_cast(n)) * one_over_x_pow_two_n_minus_one) / (n2 * (n2 - 1U)); + + if ((n >= 3U) && (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; + } + if (n > number_of_bernoullis_b2n) + return policies::raise_evaluation_error("scaled_tgamma_no_lanczos<%1%>()", "Exceeded maximum series iterations without reaching convergence, best approximation was %1%", T(exp(sum) * half_ln_two_pi_over_z), pol); + + sum += term; + + // Sanity check for divergence: + T fterm = fabs(term); + if(fterm > last_term) + return policies::raise_evaluation_error("scaled_tgamma_no_lanczos<%1%>()", "Series became divergent without reaching convergence, best approximation was %1%", T(exp(sum) * half_ln_two_pi_over_z), pol); + last_term = fterm; + } + + // Complete Stirling's approximation. + T scaled_gamma_value = islog ? T(sum + log(half_ln_two_pi_over_z)) : T(exp(sum) * half_ln_two_pi_over_z); + return scaled_gamma_value; } // Forward declaration of the lgamma_imp template specialization. @@ -391,7 +453,7 @@ T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&) // Check if the argument of tgamma is identically zero. const bool is_at_zero = (z == 0); - if((is_at_zero) || ((boost::math::isinf)(z) && (z < 0))) + if((boost::math::isnan)(z) || (is_at_zero) || ((boost::math::isinf)(z) && (z < 0))) return policies::raise_domain_error(function, "Evaluation of tgamma at %1%.", z, pol); const bool b_neg = (z < 0); @@ -422,7 +484,7 @@ T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&) // 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(); + const int min_arg_for_recursion = minimum_argument_for_bernoulli_recursion(); int n_recur; @@ -436,13 +498,20 @@ T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&) { n_recur = 0; } - - const T log_gamma_value = lgamma_imp(zz, pol, lanczos::undefined_lanczos()); - - if(log_gamma_value > tools::log_max_value()) + if (!n_recur) + { + if (zz > tools::log_max_value()) + return policies::raise_overflow_error(function, 0, pol); + if (log(zz) * zz / 2 > tools::log_max_value()) + return policies::raise_overflow_error(function, 0, pol); + } + T gamma_value = scaled_tgamma_no_lanczos(zz, pol); + T power_term = pow(zz, zz / 2); + T exp_term = exp(-zz); + gamma_value *= (power_term * exp_term); + if(!n_recur && (tools::max_value() / power_term < gamma_value)) return policies::raise_overflow_error(function, 0, pol); - - T gamma_value = exp(log_gamma_value); + gamma_value *= power_term; // Rescale the result using downward recursion if necessary. if(n_recur) @@ -500,7 +569,8 @@ inline T log_gamma_near_1(const T& z, Policy const& pol) { // // This is for the multiprecision case where there is - // no lanczos support... + // no lanczos support, use a taylor series at z = 1, + // see https://www.wolframalpha.com/input/?i=taylor+series+lgamma(x)+at+x+%3D+1 // BOOST_MATH_STD_USING // ADL of std names @@ -508,20 +578,17 @@ inline T log_gamma_near_1(const T& z, Policy const& pol) T result = -constants::euler() * z; - T power_term = z * z; - T term; - unsigned j = 0; + T power_term = z * z / 2; + int n = 2; + T term = 0; do { - term = boost::math::zeta(j + 2, pol) * power_term / (j + 2); - if(j & 1) - result -= term; - else - result += term; - power_term *= z; - ++j; - } while(fabs(result) * tools::epsilon() < fabs(term)); + term = power_term * boost::math::polygamma(n - 1, T(1)); + result += term; + ++n; + power_term *= z / n; + } while (fabs(result) * tools::epsilon() < fabs(term)); return result; } @@ -550,13 +617,15 @@ T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sig // Special case handling of small factorials: if((!b_neg) && floor_of_z_is_equal_to_z && (z < boost::math::max_factorial::value)) { + if (sign) + *sign = 1; return log(boost::math::unchecked_factorial(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(); + const int min_arg_for_recursion = minimum_argument_for_bernoulli_recursion(); T log_gamma_value; @@ -565,68 +634,41 @@ T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sig // 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 (sign) + * sign = 1; if(fabs(z - 1) < 0.25) { - return log_gamma_near_1(T(zz - 1), pol); + log_gamma_value = log_gamma_near_1(T(zz - 1), pol); } else if(fabs(z - 2) < 0.25) { - return log_gamma_near_1(T(zz - 2), pol) + log(zz - 1); + log_gamma_value = log_gamma_near_1(T(zz - 2), pol) + log(zz - 1); } else if (z > -tools::root_epsilon()) { // 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; + *sign = z < 0 ? -1 : 1; + return log(abs(gamma_imp(z, pol, lanczos::undefined_lanczos()))); } 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()))); + T g = gamma_imp(zz, pol, lanczos::undefined_lanczos()); + if (sign) + { + *sign = g < 0 ? -1 : 1; + } + log_gamma_value = log(abs(g)); } } else { // Perform the Bernoulli series expansion of Stirling's approximation. - - const std::size_t number_of_bernoullis_b2n = highest_bernoulli_index(); - - 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 * boost::math::tools::epsilon()) * 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(n * 2U); - - 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 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()) / 2; - - log_gamma_value = ((((zz - boost::math::constants::half()) * log(zz)) - zz) + half_ln_two_pi) + sum; + T sum = scaled_tgamma_no_lanczos(zz, pol, true); + log_gamma_value = zz * (log(zz) - 1) + sum; } int sign_of_result = 1; @@ -769,6 +811,8 @@ T full_igamma_prefix(T a, T z, const Policy& pol) BOOST_MATH_STD_USING T prefix; + if (z > tools::max_value()) + return 0; T alz = a * log(z); if(z >= 1) @@ -818,6 +862,8 @@ template T regularised_gamma_prefix(T a, T z, const Policy& pol, const Lanczos& l) { BOOST_MATH_STD_USING + if (z >= tools::max_value()) + return 0; T agh = a + static_cast(Lanczos::g()) - T(0.5); T prefix; T d = ((z - a) - static_cast(Lanczos::g()) + T(0.5)) / agh; @@ -897,49 +943,70 @@ T regularised_gamma_prefix(T a, T z, const Policy& pol, const Lanczos& l) // And again, without Lanczos support: // template -T regularised_gamma_prefix(T a, T z, const Policy& pol, const lanczos::undefined_lanczos&) +T regularised_gamma_prefix(T a, T z, const Policy& pol, const lanczos::undefined_lanczos& l) { BOOST_MATH_STD_USING - T limit = (std::max)(T(10), a); - T sum = detail::lower_gamma_series(a, limit, pol) / a; - sum += detail::upper_gamma_fraction(a, limit, ::boost::math::policies::get_epsilon()); - - if(a < 10) + if((a < 1) && (z < 1)) { - // special case for small a: - T prefix = pow(z / 10, a); - prefix *= exp(10-z); - if(0 == prefix) + // No overflow possible since the power terms tend to unity as a,z -> 0 + return pow(z, a) * exp(-z) / boost::math::tgamma(a, pol); + } + else if(a > minimum_argument_for_bernoulli_recursion()) + { + T scaled_gamma = scaled_tgamma_no_lanczos(a, pol); + T power_term = pow(z / a, a / 2); + T a_minus_z = a - z; + if ((0 == power_term) || (fabs(a_minus_z) > tools::log_max_value())) { - prefix = pow((z * exp((10-z)/a)) / 10, a); + // The result is probably zero, but we need to be sure: + return exp(a * log(z / a) + a_minus_z - log(scaled_gamma)); } - prefix /= sum; - return prefix; + return (power_term * exp(a_minus_z)) * (power_term / scaled_gamma); } - - T zoa = z / a; - T amz = a - z; - T alzoa = a * log(zoa); - T prefix; - if(((std::min)(alzoa, amz) <= tools::log_min_value()) || ((std::max)(alzoa, amz) >= tools::log_max_value())) + else { - T amza = amz / a; - if((amza <= tools::log_min_value()) || (amza >= tools::log_max_value())) + // + // Usual case is to calculate the prefix at a+shift and recurse down + // to the value we want: + // + const int min_z = minimum_argument_for_bernoulli_recursion(); + long shift = 1 + ltrunc(min_z - a); + T result = regularised_gamma_prefix(T(a + shift), z, pol, l); + if (result != 0) { - prefix = exp(alzoa + amz); + for (long i = 0; i < shift; ++i) + { + result /= z; + result *= a + i; + } + return result; } else { - prefix = pow(zoa * exp(amza), a); + // + // We failed, most probably we have z << 1, try again, this time + // we calculate z^a e^-z / tgamma(a+shift), combining power terms + // as we go. And again recurse down to the result. + // + T scaled_gamma = scaled_tgamma_no_lanczos(T(a + shift), pol); + T power_term_1 = pow(z / (a + shift), a); + T power_term_2 = pow(a + shift, -shift); + T power_term_3 = exp(a + shift - z); + if ((0 == power_term_1) || (0 == power_term_2) || (0 == power_term_3) || (fabs(a + shift - z) > tools::log_max_value())) + { + // We have no test case that gets here, most likely the type T + // has a high precision but low exponent range: + return exp(a * log(z) - z - boost::math::lgamma(a, pol)); + } + result = power_term_1 * power_term_2 * power_term_3 / scaled_gamma; + for (long i = 0; i < shift; ++i) + { + result *= a + i; + } + return result; } } - else - { - prefix = pow(zoa, a) * exp(amz); - } - prefix /= sum; - return prefix; } // // Upper gamma fraction for very small a: @@ -1035,6 +1102,37 @@ T finite_half_gamma_q(T a, T x, T* p_derivative, const Policy& pol) } return e; } +// +// Asymptotic approximation for large argument, see: https://dlmf.nist.gov/8.11#E2 +// +template +struct incomplete_tgamma_large_x_series +{ + typedef T result_type; + incomplete_tgamma_large_x_series(const T& a, const T& x) + : a_poch(a - 1), z(x), term(1) {} + T operator()() + { + T result = term; + term *= a_poch / z; + a_poch -= 1; + return result; + } + T a_poch, z, term; +}; + +template +T incomplete_tgamma_large_x(const T& a, const T& x, const Policy& pol) +{ + BOOST_MATH_STD_USING + incomplete_tgamma_large_x_series s(a, x); + boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations(); + T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon(), max_iter); + boost::math::policies::check_series_iterations("boost::math::tgamma<%1%>(%1%,%1%)", max_iter, pol); + return result; +} + + // // Main incomplete gamma entry point, handles all four incomplete gamma's: // @@ -1151,6 +1249,12 @@ T gamma_incomplete_imp(T a, T x, bool normalised, bool invert, { eval_method = 6; } + else if ((x > 1000) && ((a < x) || (fabs(a - 50) / x < 1))) + { + // calculate Q via asymptotic approximation: + invert = !invert; + eval_method = 7; + } else if(x < 0.5) { // @@ -1365,7 +1469,22 @@ T gamma_incomplete_imp(T a, T x, bool normalised, bool invert, // 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 (p_derivative) + *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type()); + break; } + case 7: + { + // x is large, + // Compute Q: + result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol); + if (p_derivative) + *p_derivative = result; + result /= x; + if (result != 0) + result *= incomplete_tgamma_large_x(a, x, pol); + break; + } } if(normalised && (result > 1)) @@ -1451,40 +1570,59 @@ T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const Lanczos& // And again without Lanczos support this time: // template -T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const lanczos::undefined_lanczos&) +T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const lanczos::undefined_lanczos& l) { BOOST_MATH_STD_USING + // - // 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. + // We adjust z and delta so that both z and z+delta are large enough for + // Sterling's approximation to hold. We can then calculate the ratio + // for the adjusted values, and rescale back down to z and z+delta. // - T prefix = 1; - T zd = z + delta; - while((zd < 6) && (z < 6)) + // Get the required shifts first: + // + long numerator_shift = 0; + long denominator_shift = 0; + const int min_z = minimum_argument_for_bernoulli_recursion(); + + if (min_z > z) + numerator_shift = 1 + ltrunc(min_z - z); + if (min_z > z + delta) + denominator_shift = 1 + ltrunc(min_z - z - delta); + // + // If the shifts are zero, then we can just combine scaled tgamma's + // and combine the remaining terms: + // + if (numerator_shift == 0 && denominator_shift == 0) { - prefix /= z; - prefix *= zd; - z += 1; - zd += 1; + T scaled_tgamma_num = scaled_tgamma_no_lanczos(z, pol); + T scaled_tgamma_denom = scaled_tgamma_no_lanczos(T(z + delta), pol); + T result = scaled_tgamma_num / scaled_tgamma_denom; + result *= exp(z * boost::math::log1p(-delta / (z + delta), pol)) * pow((delta + z) / constants::e(), -delta); + return result; } - if(delta < 10) + // + // We're going to have to rescale first, get the adjusted z and delta values, + // plus the ratio for the adjusted values: + // + T zz = z + numerator_shift; + T dd = delta - (numerator_shift - denominator_shift); + T ratio = tgamma_delta_ratio_imp_lanczos(zz, dd, pol, l); + // + // Use gamma recurrence relations to get back to the original + // z and z+delta: + // + for (long long i = 0; i < numerator_shift; ++i) { - prefix *= exp(-z * boost::math::log1p(delta / z, pol)); + ratio /= (z + i); + if (i < denominator_shift) + ratio *= (z + delta + i); } - else + for (long long i = numerator_shift; i < denominator_shift; ++i) { - prefix *= pow(z / zd, z); - } - prefix *= pow(constants::e() / zd, delta); - T sum = detail::lower_gamma_series(z, z, pol) / z; - sum += detail::upper_gamma_fraction(z, z, ::boost::math::policies::get_epsilon()); - T sumd = detail::lower_gamma_series(zd, zd, pol) / zd; - sumd += detail::upper_gamma_fraction(zd, zd, ::boost::math::policies::get_epsilon()); - sum /= sumd; - if(fabs(tools::max_value() / prefix) < fabs(sum)) - return policies::raise_overflow_error("boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)", "Result of tgamma is too large to represent.", pol); - return sum * prefix; + ratio *= (z + delta + i); + } + return ratio; } template