From 9086171934ca2f3ae3347752d5bf58f86f6478de Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Mon, 12 Jan 2026 18:04:17 -0800 Subject: [PATCH 1/5] Implemented large a series expansion --- .../boost/math/special_functions/gamma.hpp | 66 +++++++++++++++++++ 1 file changed, 66 insertions(+) diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp index 458dc40a3d..4353140a54 100644 --- a/include/boost/math/special_functions/gamma.hpp +++ b/include/boost/math/special_functions/gamma.hpp @@ -1290,6 +1290,33 @@ BOOST_MATH_GPU_ENABLED T incomplete_tgamma_large_x(const T& a, const T& x, const return result; } +template +struct incomplete_tgamma_lower_large_a_series +{ + typedef T result_type; + BOOST_MATH_GPU_ENABLED incomplete_tgamma_lower_large_a_series(const T& a, const T& x) + : a_poch(a + 1), z(x), term(1 / a) {} + BOOST_MATH_GPU_ENABLED T operator()() + { + T result = term; + term *= z / a_poch; + a_poch += 1; + return result; + } + T a_poch, z, term; +}; + +template +T incomplete_tgamma_lower_large_a(const T& a, const T&x, const Policy & pol) +{ + BOOST_MATH_STD_USING + incomplete_tgamma_lower_large_a_series s(a, x); + boost::math::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_p<%1%>(%1%,%1%)", max_iter, pol); + return result; +} + // // Main incomplete gamma entry point, handles all four incomplete gamma's: @@ -1813,6 +1840,45 @@ BOOST_MATH_GPU_ENABLED T lgamma_incomplete_imp(T a, T x, const Policy& pol) return log(gamma_q(a, x, pol)); } +// Calculate log of incomplete gamma function +template +BOOST_MATH_GPU_ENABLED T lgamma_incomplete_lower_imp(T a, T x, const Policy& pol) +{ + using namespace boost::math; // temporary until we're in the right namespace + + BOOST_MATH_STD_USING_CORE + + // Check for invalid inputs (a < 0 or x < 0) + 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); + + // Need to change this to be more appropriate! + if (a > 100){ + return log(detail::incomplete_tgamma_lower_large_a(a, x, pol)) + a * log(x) - x - lgamma(a, pol); + } + + // + // Can't do better than taking the log of P, but... + // + // Figure out whether we need P or Q, since if we calculate P and it's too close to unity + // we will lose 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 log(gamma_p(a, x, pol)); + return log1p(-gamma_q(a, x, pol), pol); +} + // // Ratios of two gamma functions: // From 7a8006109b5a74100a179e5f6dbae000cadf429e Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Mon, 12 Jan 2026 19:51:50 -0800 Subject: [PATCH 2/5] Added boilerplate typecasting --- .../boost/math/special_functions/gamma.hpp | 30 +++++++++++++++++-- 1 file changed, 28 insertions(+), 2 deletions(-) diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp index 4353140a54..b12ef8fb48 100644 --- a/include/boost/math/special_functions/gamma.hpp +++ b/include/boost/math/special_functions/gamma.hpp @@ -1290,6 +1290,9 @@ BOOST_MATH_GPU_ENABLED T incomplete_tgamma_large_x(const T& a, const T& x, const return result; } +// +// Asymptotic approximation for large a, see https://dlmf.nist.gov/8.11#E4 +// template struct incomplete_tgamma_lower_large_a_series { @@ -1855,8 +1858,8 @@ BOOST_MATH_GPU_ENABLED T lgamma_incomplete_lower_imp(T a, T x, const Policy& 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); - // Need to change this to be more appropriate! - if (a > 100){ + // This still needs work... + if (a > x + 50){ return log(detail::incomplete_tgamma_lower_large_a(a, x, pol)) + a * log(x) - x - lgamma(a, pol); } @@ -2520,6 +2523,29 @@ BOOST_MATH_GPU_ENABLED inline tools::promote_args_t lgamma_q(T1 a, T2 z) { return lgamma_q(a, z, policies::policy<>()); } + +template +BOOST_MATH_GPU_ENABLED inline tools::promote_args_t lgamma_p(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_lower_imp(static_cast(a), + static_cast(z), forwarding_policy()), "lgamma_p<%1%>(%1%, %1%)"); +} + +template +BOOST_MATH_GPU_ENABLED inline tools::promote_args_t lgamma_p(T1 a, T2 z) +{ + return lgamma_p(a, z, policies::policy<>()); +} // // Regularised lower incomplete gamma: // From 4e9f4410c05356f8849eaa17d25a31d68c963bd1 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Mon, 12 Jan 2026 20:01:26 -0800 Subject: [PATCH 3/5] Added lgamma_p to math_fwd.hpp --- include/boost/math/special_functions/math_fwd.hpp | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/include/boost/math/special_functions/math_fwd.hpp b/include/boost/math/special_functions/math_fwd.hpp index ca16fc8d96..6c8394ddf9 100644 --- a/include/boost/math/special_functions/math_fwd.hpp +++ b/include/boost/math/special_functions/math_fwd.hpp @@ -567,6 +567,12 @@ namespace boost template BOOST_MATH_GPU_ENABLED tools::promote_args_t lgamma_q(RT1 a, RT2 z, const Policy&); + template + BOOST_MATH_GPU_ENABLED tools::promote_args_t lgamma_p(RT1 a, RT2 z); + + template + BOOST_MATH_GPU_ENABLED tools::promote_args_t lgamma_p(RT1 a, RT2 z, const Policy&); + template BOOST_MATH_GPU_ENABLED tools::promote_args_t gamma_p(RT1 a, RT2 z); @@ -1525,6 +1531,9 @@ namespace boost \ template \ BOOST_MATH_GPU_ENABLED inline boost::math::tools::promote_args_t lgamma_q(RT1 a, RT2 z){ return boost::math::lgamma_q(a, z, Policy()); }\ +\ + template \ + BOOST_MATH_GPU_ENABLED inline boost::math::tools::promote_args_t lgamma_p(RT1 a, RT2 z){ return boost::math::lgamma_p(a, z, Policy()); }\ \ template \ BOOST_MATH_GPU_ENABLED inline boost::math::tools::promote_args_t gamma_p(RT1 a, RT2 z){ return boost::math::gamma_p(a, z, Policy()); }\ From 038cb2b8876cf06f4daa02a6742772724cab4646 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Tue, 13 Jan 2026 09:04:58 -0800 Subject: [PATCH 4/5] Deleted redundant lower gamma series --- .../boost/math/special_functions/gamma.hpp | 34 ++----------------- 1 file changed, 2 insertions(+), 32 deletions(-) diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp index b12ef8fb48..22a446a1e2 100644 --- a/include/boost/math/special_functions/gamma.hpp +++ b/include/boost/math/special_functions/gamma.hpp @@ -1290,37 +1290,6 @@ BOOST_MATH_GPU_ENABLED T incomplete_tgamma_large_x(const T& a, const T& x, const return result; } -// -// Asymptotic approximation for large a, see https://dlmf.nist.gov/8.11#E4 -// -template -struct incomplete_tgamma_lower_large_a_series -{ - typedef T result_type; - BOOST_MATH_GPU_ENABLED incomplete_tgamma_lower_large_a_series(const T& a, const T& x) - : a_poch(a + 1), z(x), term(1 / a) {} - BOOST_MATH_GPU_ENABLED T operator()() - { - T result = term; - term *= z / a_poch; - a_poch += 1; - return result; - } - T a_poch, z, term; -}; - -template -T incomplete_tgamma_lower_large_a(const T& a, const T&x, const Policy & pol) -{ - BOOST_MATH_STD_USING - incomplete_tgamma_lower_large_a_series s(a, x); - boost::math::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_p<%1%>(%1%,%1%)", max_iter, pol); - return result; -} - - // // Main incomplete gamma entry point, handles all four incomplete gamma's: // @@ -1860,7 +1829,8 @@ BOOST_MATH_GPU_ENABLED T lgamma_incomplete_lower_imp(T a, T x, const Policy& pol // This still needs work... if (a > x + 50){ - return log(detail::incomplete_tgamma_lower_large_a(a, x, pol)) + a * log(x) - x - lgamma(a, pol); + std::cout << "Using approximation" << std::endl; + return log(detail::lower_gamma_series(a, x, pol)) - log(a) + a * log(x) - x - lgamma(a, pol); } // From 9abdc80f15b71ac7ca87fd426668a3f5696c9b52 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Tue, 13 Jan 2026 09:32:20 -0800 Subject: [PATCH 5/5] Changed asymptotic selection criteria to match previous implementations --- include/boost/math/special_functions/gamma.hpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp index 22a446a1e2..267e58620d 100644 --- a/include/boost/math/special_functions/gamma.hpp +++ b/include/boost/math/special_functions/gamma.hpp @@ -1827,9 +1827,8 @@ BOOST_MATH_GPU_ENABLED T lgamma_incomplete_lower_imp(T a, T x, const Policy& 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); - // This still needs work... - if (a > x + 50){ - std::cout << "Using approximation" << std::endl; + // Taken from conditions on Line 1354, 1368, 1709 + if (((a > 4 * x) && (a >= max_factorial::value)) || (T(-0.4) / log(x) < a) || ((x * 0.75f < a) && x < T(1.1))){ return log(detail::lower_gamma_series(a, x, pol)) - log(a) + a * log(x) - x - lgamma(a, pol); }