From 24ede6fcdd73945c0d57962ca62088e99d557f72 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Wed, 3 Jul 2024 19:52:07 +0100 Subject: [PATCH] gamma.hpp coverage. --- include/boost/math/special_functions/gamma.hpp | 18 ++++++++++++------ test/test_tgamma_eh.cpp | 6 ++++++ test/test_tgamma_ratio.hpp | 16 ++++++++++++++++ 3 files changed, 34 insertions(+), 6 deletions(-) diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp index 523e1bc84..3ca4e4c63 100644 --- a/include/boost/math/special_functions/gamma.hpp +++ b/include/boost/math/special_functions/gamma.hpp @@ -1436,7 +1436,7 @@ T gamma_incomplete_imp(T a, T x, bool normalised, bool invert, optimised_invert = true; } else - init_value = 0; + init_value = 0; // LCOV_EXCL_LINE Unreachable for any "sensible" floating point type. } else init_value = 0; @@ -1546,15 +1546,18 @@ T gamma_incomplete_imp(T a, T x, bool normalised, bool invert, } if(p_derivative) { - // - // Need to convert prefix term to derivative: - // + /* + * We can never reach this: if((x < 1) && (tools::max_value() * x < *p_derivative)) { // overflow, just return an arbitrarily large value: *p_derivative = tools::max_value() / 2; } - + */ + BOOST_MATH_ASSERT((x >= 1) || (tools::max_value() * x >= *p_derivative)); + // + // Need to convert prefix term to derivative: + // *p_derivative /= x; } @@ -1800,7 +1803,10 @@ T tgamma_ratio_imp(T x, T y, const Policy& pol) // // Result will almost certainly overflow, try logs just in case: // - return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol)); + prefix = boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol); + if (prefix > boost::math::tools::log_max_value()) + return policies::raise_overflow_error("tgamma_ratio", nullptr, pol); + return exp(prefix); } // // Regular case, x and y both large and similar in magnitude: diff --git a/test/test_tgamma_eh.cpp b/test/test_tgamma_eh.cpp index d27e7e5c1..1ebb4e4ce 100644 --- a/test/test_tgamma_eh.cpp +++ b/test/test_tgamma_eh.cpp @@ -3,6 +3,8 @@ // Boost Software License, Version 1.0. (See accompanying file // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) +#define BOOST_MATH_PROMOTE_DOUBLE_POLICY false + #include "math_unit_test.hpp" #include #include @@ -13,5 +15,9 @@ int main() { CHECK_EQUAL(boost::math::tgamma(-200.5), 0.0); // triggers internal exception handling CHECK_EQUAL(boost::math::gamma_p(500.125, 1e-50), 0.0); // triggers internal exception handling + + // Lines that can only be hit when promotion to 80-bit reals is turned off + CHECK_ULP_CLOSE(boost::math::tgamma(44.0, 0.000001), 6.04152630633738356373551320685139975072645120000000000000000e52, 10); + CHECK_ULP_CLOSE(boost::math::gamma_p(1.0001, boost::math::gamma_p_inv(1.0001, 1e-200)), 1e-200, 10); return boost::math::test::report_errors(); } diff --git a/test/test_tgamma_ratio.hpp b/test/test_tgamma_ratio.hpp index 32a109290..37910223a 100644 --- a/test/test_tgamma_ratio.hpp +++ b/test/test_tgamma_ratio.hpp @@ -179,6 +179,16 @@ void test_spots(T, const char*) BOOST_CHECK_CLOSE_FRACTION(boost::math::tgamma_delta_ratio(T(ldexp(T(1), -1074)), T(200)), T(5.13282785052571536804189023927976812551830809667482691717029e-50L), tol * 10); BOOST_CHECK_CLOSE_FRACTION(boost::math::tgamma_delta_ratio(T(200), T(ldexp(T(1), -1074))), T(1), tol); } +#if LDBL_MAX_EXP >= 16384 + if(0 != ldexp(T(1), -16445)) + { + // This is denorm_min at 80 bit long double precision: + BOOST_CHECK_CLOSE_FRACTION(boost::math::tgamma_ratio(T(ldexp(T(1), -16445)), T(200)), T(6.956968705776973348795055087648739671964943862403874276047e4577L), tol * 50); + BOOST_CHECK_CLOSE_FRACTION(boost::math::tgamma_ratio(T(200), T(ldexp(T(1), -16445))), T(1.43740764446678254374798600325709882852078655978124426386e-4578L), tol * 10); + BOOST_CHECK_CLOSE_FRACTION(boost::math::tgamma_delta_ratio(T(ldexp(T(1), -16445)), T(200)), T(6.956968705776973348795055087648739671964943862403874276047e4577L), tol * 10); + BOOST_CHECK_CLOSE_FRACTION(boost::math::tgamma_delta_ratio(T(200), T(ldexp(T(1), -16445))), T(1), tol); + } +#endif } // // Coverage: @@ -186,4 +196,10 @@ void test_spots(T, const char*) BOOST_CHECK_EQUAL(boost::math::tgamma_delta_ratio(T(2), T(0)), T(1)); BOOST_CHECK_EQUAL(boost::math::tgamma_delta_ratio(T(200), T(0)), T(1)); BOOST_CHECK_EQUAL(boost::math::tgamma_delta_ratio(T(2000), T(0)), T(1)); + + BOOST_CHECK_EQUAL(boost::math::tgamma_ratio(T(0.5), T(100000)), T(0)); + BOOST_IF_CONSTEXPR(std::numeric_limits::has_infinity) + { + BOOST_CHECK_EQUAL(boost::math::tgamma_ratio(T(100000), T(0.5)), std::numeric_limits::infinity()); + } }