From d967bf73d0f7f8ace5011a09c6dd2650dd3be348 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Wed, 30 Aug 2023 11:50:28 +0100 Subject: [PATCH] Improve accuracy in ibeta power terms usig Sterling's approximation (the non Lanczos case). Completes work started here: https://github.com/boostorg/math/pull/1007 --- .../detail/inv_discrete_quantile.hpp | 2 +- include/boost/math/special_functions/beta.hpp | 102 ++++++++++++++++-- test/test_binomial.cpp | 5 +- 3 files changed, 97 insertions(+), 12 deletions(-) diff --git a/include/boost/math/distributions/detail/inv_discrete_quantile.hpp b/include/boost/math/distributions/detail/inv_discrete_quantile.hpp index 7ded3d39ad..2f5cae7432 100644 --- a/include/boost/math/distributions/detail/inv_discrete_quantile.hpp +++ b/include/boost/math/distributions/detail/inv_discrete_quantile.hpp @@ -147,7 +147,7 @@ typename Dist::value_type // we're assuming that "guess" is likely to be accurate // to the nearest int or so: // - else if(adder != 0) + else if((adder != 0) && (a + adder != a)) { // // If we're looking for a large result, then bump "adder" up diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp index 8121c1622d..73b9e3d52c 100644 --- a/include/boost/math/special_functions/beta.hpp +++ b/include/boost/math/special_functions/beta.hpp @@ -473,20 +473,108 @@ T ibeta_power_terms(T a, if ((shift_a == 0) && (shift_b == 0)) { T power1, power2; + bool need_logs = false; if (a < b) { - power1 = pow((x * y * c * c) / (a * b), a); - power2 = pow((y * c) / b, b - a); + BOOST_IF_CONSTEXPR(std::numeric_limits::has_infinity) + { + power1 = pow((x * y * c * c) / (a * b), a); + power2 = pow((y * c) / b, b - a); + } + else + { + // We calculate these logs purely so we can check for overflow in the power functions + T l1 = log((x * y * c * c) / (a * b)); + T l2 = log((y * c) / b); + if ((l1 * a > tools::log_min_value()) && (l1 * a < tools::log_max_value()) && (l2 * (b - a) < tools::log_max_value()) && (l2 * (b - a) > tools::log_min_value())) + { + power1 = pow((x * y * c * c) / (a * b), a); + power2 = pow((y * c) / b, b - a); + } + else + { + need_logs = true; + } + } } else { - power1 = pow((x * y * c * c) / (a * b), b); - power2 = pow((x * c) / a, a - b); + BOOST_IF_CONSTEXPR(std::numeric_limits::has_infinity) + { + power1 = pow((x * y * c * c) / (a * b), b); + power2 = pow((x * c) / a, a - b); + } + else + { + // We calculate these logs purely so we can check for overflow in the power functions + T l1 = log((x * y * c * c) / (a * b)) * b; + T l2 = log((x * c) / a) * (a - b); + if ((l1 * a > tools::log_min_value()) && (l1 * a < tools::log_max_value()) && (l2 * (b - a) < tools::log_max_value()) && (l2 * (b - a) > tools::log_min_value())) + { + power1 = pow((x * y * c * c) / (a * b), b); + power2 = pow((x * c) / a, a - b); + } + else + need_logs = true; + } } - if (!(boost::math::isnormal)(power1) || !(boost::math::isnormal)(power2)) + BOOST_IF_CONSTEXPR(std::numeric_limits::has_infinity) { - // We have to use logs :( - return prefix * exp(a * log(x * c / a) + b * log(y * c / b)) * scaled_tgamma_no_lanczos(c, pol) / (scaled_tgamma_no_lanczos(a, pol) * scaled_tgamma_no_lanczos(b, pol)); + if (!(boost::math::isnormal)(power1) || !(boost::math::isnormal)(power2)) + { + need_logs = true; + } + } + if (need_logs) + { + // + // We want: + // + // (xc / a)^a (yc / b)^b + // + // But we know that one or other term will over / underflow and combining the logs will be next to useless as that will cause significant cancellation. + // If we assume b > a and express z ^ b as(z ^ b / a) ^ a with z = (yc / b) then we can move one power term inside the other : + // + // ((xc / a) * (yc / b)^(b / a))^a + // + // However, we're not quite there yet, as the term being exponentiated is quite likely to be close to unity, so let: + // + // xc / a = 1 + (xb - ya) / a + // + // analogously let : + // + // 1 + p = (yc / b) ^ (b / a) = 1 + expm1((b / a) * log1p((ya - xb) / b)) + // + // so putting the two together we have : + // + // exp(a * log1p((xb - ya) / a + p + p(xb - ya) / a)) + // + // Analogously, when a > b we can just swap all the terms around. + // + // Finally, there are a few cases (x or y is unity) when the above logic can't be used + // or where there is no logarithmic cancellation and accuracy is better just using + // the regular formula: + // + T xc_a = x * c / a; + T yc_b = y * c / b; + if ((x == 1) || (y == 1) || ((fabs(xc_a - 1) > 0.25) && (fabs(yc_b - 1) > 0.25))) + { + // The above logic fails, the result is almost certainly zero: + power1 = exp(log(xc_a) * a + log(yc_b) * b); + power2 = 1; + } + else if (b > a) + { + T p = boost::math::expm1((b / a) * boost::math::log1p((y * a - x * b) / b)); + power1 = exp(a * log1p((x * b - y * a) / a + p * (x * c / a))); + power2 = 1; + } + else + { + T p = boost::math::expm1((a / b) * boost::math::log1p((x * b - y * a) / a)); + power1 = exp(b * log1p((y * a - x * b) / b + p * (y * c / b))); + power2 = 1; + } } return prefix * power1 * power2 * scaled_tgamma_no_lanczos(c, pol) / (scaled_tgamma_no_lanczos(a, pol) * scaled_tgamma_no_lanczos(b, pol)); } diff --git a/test/test_binomial.cpp b/test/test_binomial.cpp index 34baeaa6be..6a8f077f34 100644 --- a/test/test_binomial.cpp +++ b/test/test_binomial.cpp @@ -716,10 +716,7 @@ void test_spots(RealType T) check_out_of_range >(1, 1); // (All) valid constructor parameter values. - // TODO: Generic ibeta_power_terms has accuracy issue when long - // double is not precise enough, causing overflow in this case. - if(!std::is_same::value || - sizeof(boost::math::concepts::real_concept_base_type) > 8) + // See bug reported here: https://github.com/boostorg/math/pull/1007 { using namespace boost::math::policies; typedef policy > Policy;