Skip to content

Keep lgamma_p finite where P underflows, and add ulps plots to the Python layer - #1565

Merged
NAThompson merged 1 commit into
developfrom
lgamma-p-underflow
Oct 12, 2026
Merged

NAThompson merged 1 commit into
developfrom
lgamma-p-underflow

Conversation

@NAThompson

Copy link
Copy Markdown
Collaborator

lgamma_p works in logs (a log-space Lanczos prefix) for large a and wherever P underflows; entropy for the arcsine and Bernoulli distributions; boost.math.ulps draws 1D and 2D ulps plots against mpmath.

Closes #1564

lgamma_p

lgamma_p used to take log(gamma_p(a, x)). Once P underflowed, that gave log(0) = -inf, which was reported as an overflow, and log P lost accuracy as P went subnormal. Now, for x < a, it uses method 2 in logs, log P = log(prefix) + log(series / a), in two cases:

  • P underflows: for every type.
  • a > 1000: for types with a Lanczos approximation, because log P is much better conditioned than P there.

The log of the prefix x^a e^-x / Γ(a) is a·log1pmx(d) + x(½ − g)/agh + ½·log(agh/e) − log(lanczos_sum_expG_scaled(a)), the logarithm of regularised_gamma_prefix's large-a form. It has none of the cancellation in a·log(x) − x − lgamma(a), which had reached 540 eps at a = 10⁷. Elsewhere log(gamma_p) stays, since it is the more accurate of the two for smaller a.

Double, 300 × 300 grid, against mpmath (x < a) and Arb (x ≥ a):

develop

this PR

  • Spurious overflows: the cyan region (8102 points) is gone, along with the 10¹² ulps band where P is subnormal.
  • Below x = a: for x < a the worst case is now about 8 ulps.
  • Not changed here: the bright band just above x = a for large a is log1p(-gamma_q(a, x)). There Q's relative error reaches 10⁶ eps (e.g. a = 3.8·10⁶, x = 1.016a), which is a gamma_q accuracy issue unchanged from develop.

Also

  • entropy(arcsine_distribution) = log(π(x_max − x_min)/4).
  • entropy(bernoulli_distribution) = −p log p − (1 − p) log1p(−p).
  • boost.math.ulps.plot / heatmap (nanobind layer) classify each point as ulps, spurious overflow (cyan), NaN or domain error (green), or out of range and handled (gray). The plots above were made with it.

🤖 Generated with Claude Code

…thon layer

lgamma_p works in logs (a log-space Lanczos prefix) for large a and wherever P underflows; entropy for the arcsine and Bernoulli distributions; boost.math.ulps draws 1D and 2D ulps plots against mpmath.
Closes #1564

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@NAThompson
NAThompson merged commit 98194df into develop Oct 12, 2026
45 of 50 checks passed
@NAThompson
NAThompson deleted the lgamma-p-underflow branch October 12, 2026 05:20
@ckormanyos

ckormanyos commented Oct 12, 2026 •

Copy link
Copy Markdown
Member

Hi Nick (@NAThompson) I don't usually flag too many warnings with Math, but I picked up some -Wshadow in one of my other projects. I'm not entirely sure, but I think it was related to this PR (or a similar one).

Could you eliminate the Wshadow warnings please? I think shadowing is one of those things clients will want to avoid.

In file included from ../boost-root/boost/multiprecision/cpp_bin_float.hpp:36,
                 from ./test/independent_algebra_test_decwide_t_boost_cpp.h:47,
                 from ./test/independent_algebra_test_decwide_t.h:13,
                 from test/test_decwide_t_algebra.cpp:61:
../boost-root/boost/math/special_functions/gamma.hpp: In function ‘T boost::math::detail::regularised_gamma_prefix(T, T, const Policy&, const boost::math::lanczos::undefined_lanczos&)’:
../boost-root/boost/math/special_functions/gamma.hpp:1146:12: error: declaration of ‘T l’ shadows a parameter [-Werror=shadow]
 1146 |          T l = a * BOOST_MATH_NAMESPACE::log1pmx(d, pol);
      |            ^
../boost-root/boost/math/special_functions/gamma.hpp:1129:91: note: shadowed declaration is here
 1129 | T regularised_gamma_prefix(T a, T z, const Policy& pol, const lanczos::undefined_lanczos& l)
      |                                                         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^
cc1plus: all warnings being treated as errors
In file included from examples/example009a_boost_math_standalone.cpp:36:
../boost-root/boost/math/special_functions/gamma.hpp: In function ‘T boost::math::detail::regularised_gamma_prefix(T, T, const Policy&, const boost::math::lanczos::undefined_lanczos&)’:
../boost-root/boost/math/special_functions/gamma.hpp:1146:12: error: declaration of ‘T l’ shadows a parameter [-Werror=shadow]
 1146 |          T l = a * BOOST_MATH_NAMESPACE::log1pmx(d, pol);
      |            ^
../boost-root/boost/math/special_functions/gamma.hpp:1129:91: note: shadowed declaration is here
 1129 | T regularised_gamma_prefix(T a, T z, const Policy& pol, const lanczos::undefined_lanczos& l)
      |                                                         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^
cc1plus: all warnings being treated as errors
In file included from examples/example009b_boost_math_standalone.cpp:45:
../boost-root/boost/math/special_functions/gamma.hpp: In function ‘T boost::math::detail::regularised_gamma_prefix(T, T, const Policy&, const boost::math::lanczos::undefined_lanczos&)’:
../boost-root/boost/math/special_functions/gamma.hpp:1146:12: error: declaration of ‘T l’ shadows a parameter [-Werror=shadow]
 1146 |          T l = a * BOOST_MATH_NAMESPACE::log1pmx(d, pol);
      |            ^
../boost-root/boost/math/special_functions/gamma.hpp:1129:91: note: shadowed declaration is here
 1129 | T regularised_gamma_prefix(T a, T z, const Policy& pol, const lanczos::undefined_lanczos& l)
      |                                                         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^
cc1plus: all warnings being treated as errors
make: [Makefile:123: compile_and_link] Error 1 (ignored)

@NAThompson

Copy link
Copy Markdown
Collaborator Author

@ckormanyos : Fixing in the entropy PR to preserve CI.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

lgamma_p throws overflow when P underflows

2 participants