From d4afedc638cb34690365cbbb8853545c4ad1e77f Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Sun, 28 Dec 2025 14:06:20 -0800 Subject: [PATCH 1/3] boost::math::detail log incomplete gamma function implemented --- .../boost/math/special_functions/gamma.hpp | 45 +++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp index 47b4f3d68d..037c7cba02 100644 --- a/include/boost/math/special_functions/gamma.hpp +++ b/include/boost/math/special_functions/gamma.hpp @@ -1772,6 +1772,51 @@ BOOST_MATH_GPU_ENABLED T gamma_incomplete_imp(T a, T x, bool normalised, bool in return gamma_incomplete_imp_final(T(a), T(x), normalised, invert, pol, p_derivative); } +// Calculate log of incomplete gamma function +template +T lgamma_incomplete_imp_final(T a, T x, const Policy& pol) +{ + using namespace boost::math; // temporary until we're in the right namespace + + BOOST_MATH_STD_USING_CORE + + if (((x > 1000) && ((a < x) || (fabs(a - 50) / x < 1))) || ((x > tools::log_max_value() - 10) && (x > a))) + { + // + // Take the logarithmic version of the asymtotic expansion: + // + return log(detail::incomplete_tgamma_large_x(a, x, pol)) + a * log(x) - x - lgamma(a, pol) - log(x); + } + // + // Can't do better than taking the log of Q, but... + // + // Figure out whether we need P or Q, since if we calculate Q and it's too close to unity + // we will loose precision in the result, selection logic here is extracted from gamma_incomplete_imp_final: + // + bool need_p = false; + if ((x < 0.5) && (T(-0.4) / log(x) < a)) + need_p = true; + else if ((x < 1.1) && (x >= 0.5) && (x * 0.75f < a)) + need_p = true; + else if ((x < a) && (x >= 1.1)) + need_p = true; + + if (need_p) + return log1p(-gamma_p(a, x, pol), pol); + return log(gamma_q(a, x, pol)); +} + +template +T lgamma_incomplete_imp(T a, T x, const Policy& pol){ + constexpr auto function = "boost::math::lgamma_p<%1%>(%1%, %1%)"; + if(a <= 0) + return policies::raise_domain_error(function, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol); + if(x < 0) + return policies::raise_domain_error(function, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol); + + // If input is valid proceed as normal + return lgamma_incomplete_imp_final(T(a), T(x), pol) +} // // Ratios of two gamma functions: // From 49950ef147437e660637be0e4b28f224f4c8a4c9 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Sun, 28 Dec 2025 14:47:50 -0800 Subject: [PATCH 2/3] Implemented log incomplete gamma with/without policy [skip ci] --- .../boost/math/special_functions/gamma.hpp | 25 ++++++++++++++++++- 1 file changed, 24 insertions(+), 1 deletion(-) diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp index 037c7cba02..a8347494b5 100644 --- a/include/boost/math/special_functions/gamma.hpp +++ b/include/boost/math/special_functions/gamma.hpp @@ -1808,7 +1808,7 @@ T lgamma_incomplete_imp_final(T a, T x, const Policy& pol) template T lgamma_incomplete_imp(T a, T x, const Policy& pol){ - constexpr auto function = "boost::math::lgamma_p<%1%>(%1%, %1%)"; + constexpr auto function = "boost::math::lgamma_q<%1%>(%1%, %1%)"; if(a <= 0) return policies::raise_domain_error(function, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol); if(x < 0) @@ -2435,6 +2435,29 @@ BOOST_MATH_GPU_ENABLED inline tools::promote_args_t { return gamma_q(a, z, policies::policy<>()); } + +template +inline tools::promote_args_t lgamma_q(T1 a, T2 z, const Policy& /* pol */) +{ + typedef tools::promote_args_t result_type; + typedef typename policies::evaluation::type value_type; + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + + return policies::checked_narrowing_cast( + detail::lgamma_incomplete_imp(static_cast(a), + static_cast(z), forwarding_policy()), "lgamma_q<%1%>(%1%, %1%)"); +} + +template +inline tools::promote_args_t lgamma_q(T1 a, T2 z) +{ + return lgamma_q(a, z, policies::policy<>()); +} // // Regularised lower incomplete gamma: // From 4cdc487bf7999681fa0db7943e92e990c7da15af Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Sun, 28 Dec 2025 16:57:55 -0800 Subject: [PATCH 3/3] Small bug fixes [skip ci] --- include/boost/math/special_functions/gamma.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp index a8347494b5..cb25118406 100644 --- a/include/boost/math/special_functions/gamma.hpp +++ b/include/boost/math/special_functions/gamma.hpp @@ -1815,7 +1815,7 @@ T lgamma_incomplete_imp(T a, T x, const Policy& pol){ return policies::raise_domain_error(function, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol); // If input is valid proceed as normal - return lgamma_incomplete_imp_final(T(a), T(x), pol) + return lgamma_incomplete_imp_final(T(a), T(x), pol); } // // Ratios of two gamma functions: @@ -2453,7 +2453,7 @@ inline tools::promote_args_t lgamma_q(T1 a, T2 z, const Policy& /* pol * static_cast(z), forwarding_policy()), "lgamma_q<%1%>(%1%, %1%)"); } -template +template inline tools::promote_args_t lgamma_q(T1 a, T2 z) { return lgamma_q(a, z, policies::policy<>());