From 9cdc9847a64a5424f695c241cf38c9419852b775 Mon Sep 17 00:00:00 2001 From: rahulb0802 Date: Tue, 25 Aug 2026 21:49:29 -0500 Subject: [PATCH 1/4] Ports Loader (2000) SP approx to prevent cancellation in extreme parameter regimes --- include/boost/math/distributions/poisson.hpp | 48 ++++++++++++++- test/test_poisson.cpp | 61 +++++++++++++++++++- 2 files changed, 107 insertions(+), 2 deletions(-) diff --git a/include/boost/math/distributions/poisson.hpp b/include/boost/math/distributions/poisson.hpp index c5e4404335..e8a31c6b62 100644 --- a/include/boost/math/distributions/poisson.hpp +++ b/include/boost/math/distributions/poisson.hpp @@ -50,6 +50,7 @@ #include // factorials. #include // for root finding. #include +#include namespace boost { @@ -141,6 +142,51 @@ namespace boost return true; } // bool check_dist_and_prob + template + BOOST_MATH_GPU_ENABLED inline RealType stirlerr(const RealType& n) { + BOOST_MATH_STD_USING // for ADL of std functions. + using boost::math::lgamma; + + const RealType S0 = RealType(1)/12; + const RealType S1 = RealType(1)/360; + const RealType S2 = RealType(1)/1260; + const RealType S3 = RealType(1)/1680; + const RealType S4 = RealType(1)/1188; + + bool is_small = n < 15; + if (is_small) { + return lgamma(n + 1) - (n * log(n) - n + 0.5 * log(2 * boost::math::constants::pi() * n)); + } else { + RealType n2 = n * n; + return (S0 - (S1 - (S2 - (S3 - S4/n2)/n2)/n2)/n2)/n; + } + + } + + template + BOOST_MATH_GPU_ENABLED inline RealType bd0(const RealType& mean, const RealType& k) { + BOOST_MATH_STD_USING // for ADL of std functions. + + bool is_close = abs(k - mean) < RealType(0.1) * (k + mean); + + if (is_close) { + RealType v = (k - mean) / (k + mean); + RealType v2 = v * v; + RealType series_term = ((k - mean) * (k - mean)) / (k + mean); + + RealType term = 2 * k * v; + for (int i = 1; i < 11; ++i) { + term *= v2; + series_term += term / (2 * i + 1); + } + return series_term; + } else { + RealType direct = (k == 0) ? RealType(0) : k * log(k / mean) + mean - k; + return direct; + } + + } + } // namespace poisson_detail BOOST_MATH_EXPORT template > @@ -304,7 +350,7 @@ namespace boost // Special case where k and lambda are both positive if(k > 0 && mean > 0) { - return -lgamma(k+1) + k*log(mean) - mean; + return -poisson_detail::stirlerr(k) - poisson_detail::bd0(mean, k) - RealType(0.5) * log(2 * boost::math::constants::pi() * k); } result = log(pdf(dist, k)); diff --git a/test/test_poisson.cpp b/test/test_poisson.cpp index 96e5f12d73..8a55028f21 100644 --- a/test/test_poisson.cpp +++ b/test/test_poisson.cpp @@ -244,7 +244,66 @@ void test_spots(RealType) static_cast(20)), // K>> mean log(static_cast(8.277463646553730E-009)), // probability. tolerance); - + + BOOST_CHECK_CLOSE( + logpdf(poisson_distribution(static_cast(14)), // mean 14. + static_cast(14)), + static_cast(-2.244418568125061), // probability. + tolerance); + + BOOST_CHECK_CLOSE( + logpdf(poisson_distribution(static_cast(20)), // mean 20. + static_cast(18)), + static_cast(-2.472264284061216), // probability. + tolerance); + + // Cases below require around 15+ significant decimal digits to represent + // k / mean meaningfully, so skip for float. + if (std::numeric_limits::digits10 > 15) + { + BOOST_CHECK_CLOSE( + logpdf(poisson_distribution(static_cast(1000000)), + static_cast(1300000)), + static_cast(-41081.501683746894), + tolerance); + + BOOST_CHECK_CLOSE( + logpdf(poisson_distribution(static_cast(1e8)), + static_cast(8e7)), + static_cast(-2148525.91257035), + tolerance); + + BOOST_CHECK_CLOSE( + logpdf(poisson_distribution(static_cast(1e10)), + static_cast(105e9)), + static_cast(-151894402015.7727), + tolerance); + + BOOST_CHECK_CLOSE( + logpdf(poisson_distribution(static_cast(1e15)), + static_cast(8e14)), // |v| > 0.1 boundary + static_cast(-21485158948650.273), + tolerance); + + BOOST_CHECK_CLOSE( + logpdf(poisson_distribution(static_cast(1e15)), + static_cast(12e14)), // |v| < 0.1 boundary + static_cast(-18785868152763.832), + tolerance); + + BOOST_CHECK_CLOSE( + logpdf(poisson_distribution(static_cast(1e16)), + static_cast(1e16)), // old formula returns 0.0 here + static_cast(-19.339619277157038), + tolerance); + + BOOST_CHECK_CLOSE( + logpdf(poisson_distribution(static_cast(5e15)), + static_cast(5e15)), + static_cast(-18.993045686877064), + tolerance); + } + // CDF BOOST_CHECK_CLOSE( cdf(poisson_distribution(static_cast(1)), // mean unity. From 82a1da8785b4986be36dba00f0d07bd07d01d7c5 Mon Sep 17 00:00:00 2001 From: rahulb0802 Date: Tue, 25 Aug 2026 21:55:36 -0500 Subject: [PATCH 2/4] Add comments to approximation code and testing --- include/boost/math/distributions/poisson.hpp | 7 ++++++- test/test_poisson.cpp | 18 +++++++++--------- 2 files changed, 15 insertions(+), 10 deletions(-) diff --git a/include/boost/math/distributions/poisson.hpp b/include/boost/math/distributions/poisson.hpp index e8a31c6b62..68c0463eb6 100644 --- a/include/boost/math/distributions/poisson.hpp +++ b/include/boost/math/distributions/poisson.hpp @@ -147,12 +147,14 @@ namespace boost BOOST_MATH_STD_USING // for ADL of std functions. using boost::math::lgamma; + // Stirling's series coefficients const RealType S0 = RealType(1)/12; const RealType S1 = RealType(1)/360; const RealType S2 = RealType(1)/1260; const RealType S3 = RealType(1)/1680; const RealType S4 = RealType(1)/1188; + // Use Stirling's series if n is small; use the direct formula otherwise bool is_small = n < 15; if (is_small) { return lgamma(n + 1) - (n * log(n) - n + 0.5 * log(2 * boost::math::constants::pi() * n)); @@ -167,9 +169,10 @@ namespace boost BOOST_MATH_GPU_ENABLED inline RealType bd0(const RealType& mean, const RealType& k) { BOOST_MATH_STD_USING // for ADL of std functions. + // Calculate v = (k - mean) / (k + mean) from Loader (2000) approximation bool is_close = abs(k - mean) < RealType(0.1) * (k + mean); - if (is_close) { + if (is_close) { // Use the series approximation if |v| < 0.1 RealType v = (k - mean) / (k + mean); RealType v2 = v * v; RealType series_term = ((k - mean) * (k - mean)) / (k + mean); @@ -181,6 +184,7 @@ namespace boost } return series_term; } else { + // Use the direct formula if |v| >= 0.1 RealType direct = (k == 0) ? RealType(0) : k * log(k / mean) + mean - k; return direct; } @@ -350,6 +354,7 @@ namespace boost // Special case where k and lambda are both positive if(k > 0 && mean > 0) { + // Use the Loader (2000) saddle-point approximation for logpdf calculation return -poisson_detail::stirlerr(k) - poisson_detail::bd0(mean, k) - RealType(0.5) * log(2 * boost::math::constants::pi() * k); } diff --git a/test/test_poisson.cpp b/test/test_poisson.cpp index 8a55028f21..9147e6787f 100644 --- a/test/test_poisson.cpp +++ b/test/test_poisson.cpp @@ -248,13 +248,13 @@ void test_spots(RealType) BOOST_CHECK_CLOSE( logpdf(poisson_distribution(static_cast(14)), // mean 14. static_cast(14)), - static_cast(-2.244418568125061), // probability. + static_cast(-2.244418568125061), // probability (already in log space). tolerance); BOOST_CHECK_CLOSE( logpdf(poisson_distribution(static_cast(20)), // mean 20. static_cast(18)), - static_cast(-2.472264284061216), // probability. + static_cast(-2.472264284061216), // probability (already in log space). tolerance); // Cases below require around 15+ significant decimal digits to represent @@ -262,43 +262,43 @@ void test_spots(RealType) if (std::numeric_limits::digits10 > 15) { BOOST_CHECK_CLOSE( - logpdf(poisson_distribution(static_cast(1000000)), + logpdf(poisson_distribution(static_cast(1000000)), // mean 1000000. static_cast(1300000)), static_cast(-41081.501683746894), tolerance); BOOST_CHECK_CLOSE( - logpdf(poisson_distribution(static_cast(1e8)), + logpdf(poisson_distribution(static_cast(1e8)), // mean 1e8. static_cast(8e7)), static_cast(-2148525.91257035), tolerance); BOOST_CHECK_CLOSE( - logpdf(poisson_distribution(static_cast(1e10)), + logpdf(poisson_distribution(static_cast(1e10)), // mean 1e10. static_cast(105e9)), static_cast(-151894402015.7727), tolerance); BOOST_CHECK_CLOSE( - logpdf(poisson_distribution(static_cast(1e15)), + logpdf(poisson_distribution(static_cast(1e15)), // mean 1e15. static_cast(8e14)), // |v| > 0.1 boundary static_cast(-21485158948650.273), tolerance); BOOST_CHECK_CLOSE( - logpdf(poisson_distribution(static_cast(1e15)), + logpdf(poisson_distribution(static_cast(1e15)), // mean 1e15. static_cast(12e14)), // |v| < 0.1 boundary static_cast(-18785868152763.832), tolerance); BOOST_CHECK_CLOSE( - logpdf(poisson_distribution(static_cast(1e16)), + logpdf(poisson_distribution(static_cast(1e16)), // mean 1e16. static_cast(1e16)), // old formula returns 0.0 here static_cast(-19.339619277157038), tolerance); BOOST_CHECK_CLOSE( - logpdf(poisson_distribution(static_cast(5e15)), + logpdf(poisson_distribution(static_cast(5e15)), // mean 5e15. static_cast(5e15)), static_cast(-18.993045686877064), tolerance); From 487009bc32d72892c68fe1331905f2351973cde7 Mon Sep 17 00:00:00 2001 From: rahulb0802 Date: Tue, 25 Aug 2026 22:06:51 -0500 Subject: [PATCH 3/4] mpmath reference in testing comment --- test/test_poisson.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/test_poisson.cpp b/test/test_poisson.cpp index 9147e6787f..22daed882f 100644 --- a/test/test_poisson.cpp +++ b/test/test_poisson.cpp @@ -245,6 +245,8 @@ void test_spots(RealType) log(static_cast(8.277463646553730E-009)), // probability. tolerance); + // New test cases for Loader (2000) saddle-point approximation. Probs already + // in log space. Values calculated using mpmath (1000-digit precision). BOOST_CHECK_CLOSE( logpdf(poisson_distribution(static_cast(14)), // mean 14. static_cast(14)), From 369fd7dc4dc1da20b28c868fc7ac097b74fcbd2d Mon Sep 17 00:00:00 2001 From: rahulb0802 Date: Fri, 28 Aug 2026 10:19:59 -0500 Subject: [PATCH 4/4] Use adaptive series for stirlerr and bd0 instead of fixed truncation Now we extract the Bernoulli series loop from scaled_tgamma_no_lanczos into a new helper which is reused by both functions. This is used instead of the fixed 5-term series, and bd0 also uses a convergence check rather than a fixed number of 10 terms. Also fixes a bug where k = mean caused convergence check to never trigger. We already know series is trivially zero so just skip the loop. --- include/boost/math/distributions/poisson.hpp | 51 +++++++++++-------- .../boost/math/special_functions/gamma.hpp | 29 +++++++---- 2 files changed, 49 insertions(+), 31 deletions(-) diff --git a/include/boost/math/distributions/poisson.hpp b/include/boost/math/distributions/poisson.hpp index 68c0463eb6..792c4a6fd2 100644 --- a/include/boost/math/distributions/poisson.hpp +++ b/include/boost/math/distributions/poisson.hpp @@ -142,31 +142,23 @@ namespace boost return true; } // bool check_dist_and_prob - template - BOOST_MATH_GPU_ENABLED inline RealType stirlerr(const RealType& n) { + template + BOOST_MATH_GPU_ENABLED inline RealType stirlerr(const RealType& n, const Policy& pol) { BOOST_MATH_STD_USING // for ADL of std functions. using boost::math::lgamma; - // Stirling's series coefficients - const RealType S0 = RealType(1)/12; - const RealType S1 = RealType(1)/360; - const RealType S2 = RealType(1)/1260; - const RealType S3 = RealType(1)/1680; - const RealType S4 = RealType(1)/1188; - - // Use Stirling's series if n is small; use the direct formula otherwise - bool is_small = n < 15; - if (is_small) { - return lgamma(n + 1) - (n * log(n) - n + 0.5 * log(2 * boost::math::constants::pi() * n)); - } else { - RealType n2 = n * n; - return (S0 - (S1 - (S2 - (S3 - S4/n2)/n2)/n2)/n2)/n; + // Use the direct formula for small n + if (n < boost::math::detail::minimum_argument_for_bernoulli_recursion()) { + return lgamma(n + 1) - (n * log(n) - n + RealType(0.5) * log(RealType(2) * boost::math::constants::pi() * n)); } + // Use the Stirling series approximation for large n + return boost::math::detail::bernoulli_stirling_series(n, pol); + } - template - BOOST_MATH_GPU_ENABLED inline RealType bd0(const RealType& mean, const RealType& k) { + template + BOOST_MATH_GPU_ENABLED inline RealType bd0(const RealType& mean, const RealType& k, const Policy& pol) { BOOST_MATH_STD_USING // for ADL of std functions. // Calculate v = (k - mean) / (k + mean) from Loader (2000) approximation @@ -178,9 +170,24 @@ namespace boost RealType series_term = ((k - mean) * (k - mean)) / (k + mean); RealType term = 2 * k * v; - for (int i = 1; i < 11; ++i) { - term *= v2; - series_term += term / (2 * i + 1); + + // Series is trivially zero when k = mean, so skip the loop + if (term != RealType(0)) { + RealType target_epsilon = abs(term) * boost::math::tools::epsilon(); + const boost::math::size_t max_iterations = policies::get_max_series_iterations(); + + for (boost::math::size_t i = 1U;; ++i) { + term *= v2; + RealType next = term / (2 * i + 1); + series_term += next; + // Break if the next term is less than the target epsilon + if (abs(next) < target_epsilon) { + break; + } + if (i > max_iterations) { + return policies::raise_evaluation_error("bd0<%1%>()", "Series did not converge in the allotted iterations, best approximation was %1%", series_term, pol); + } + } } return series_term; } else { @@ -355,7 +362,7 @@ namespace boost if(k > 0 && mean > 0) { // Use the Loader (2000) saddle-point approximation for logpdf calculation - return -poisson_detail::stirlerr(k) - poisson_detail::bd0(mean, k) - RealType(0.5) * log(2 * boost::math::constants::pi() * k); + return -poisson_detail::stirlerr(k, Policy()) - poisson_detail::bd0(mean, k, Policy()) - RealType(0.5) * log(RealType(2) * boost::math::constants::pi() * k); } result = log(pdf(dist, k)); diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp index d549158791..6016a4d77b 100644 --- a/include/boost/math/special_functions/gamma.hpp +++ b/include/boost/math/special_functions/gamma.hpp @@ -487,14 +487,9 @@ int minimum_argument_for_bernoulli_recursion() } template -T scaled_tgamma_no_lanczos(const T& z, const Policy& pol, bool islog = false) -{ +T bernoulli_stirling_series(const T& z, const Policy& pol) { 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_MATH_ASSERT(minimum_argument_for_bernoulli_recursion() <= z); // Perform the Bernoulli series expansion of Stirling's approximation. @@ -527,7 +522,7 @@ T scaled_tgamma_no_lanczos(const T& z, const Policy& pol, bool islog = false) } if (n > number_of_bernoullis_b2n) // Safety net, we hope to never get here: - 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); // LCOV_EXCL_LINE + return policies::raise_evaluation_error("bernoulli_stirling_series<%1%>()", "Exceeded maximum series iterations without reaching convergence, best approximation was %1%", sum, pol); // LCOV_EXCL_LINE sum += term; @@ -535,9 +530,25 @@ T scaled_tgamma_no_lanczos(const T& z, const Policy& pol, bool islog = false) T fterm = fabs(term); if(fterm > last_term) // Safety net, we hope to never get here: - 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); // LCOV_EXCL_LINE + return policies::raise_evaluation_error("bernoulli_stirling_series<%1%>()", "Series became divergent without reaching convergence, best approximation was %1%", sum, pol); // LCOV_EXCL_LINE last_term = fterm; } + return sum; +} + +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_MATH_ASSERT(minimum_argument_for_bernoulli_recursion() <= z); + + T sum = bernoulli_stirling_series(z, pol); + const T half_ln_two_pi_over_z = sqrt(boost::math::constants::two_pi() / z); // 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);