From b22b191fc8dc924a23300d589d0727d9f87e2dbb Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Sun, 11 Oct 2026 22:17:10 -0700 Subject: [PATCH] Keep lgamma_p finite where P underflows, and add ulps plots to the Python 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 --- include/boost/math/distributions/arcsine.hpp | 21 +++ .../boost/math/distributions/bernoulli.hpp | 12 ++ .../boost/math/special_functions/gamma.hpp | 41 ++++- nanobind/CMakeLists.txt | 6 +- nanobind/boost/math/ulps.py | 151 ++++++++++++++++++ nanobind/tests/test_ulps.py | 18 +++ test/test_arcsine.cpp | 4 + test/test_bernoulli.cpp | 10 ++ test/test_igamma.hpp | 12 ++ 9 files changed, 272 insertions(+), 3 deletions(-) create mode 100644 nanobind/boost/math/ulps.py create mode 100644 nanobind/tests/test_ulps.py diff --git a/include/boost/math/distributions/arcsine.hpp b/include/boost/math/distributions/arcsine.hpp index e5f286fea0..d476002fb9 100644 --- a/include/boost/math/distributions/arcsine.hpp +++ b/include/boost/math/distributions/arcsine.hpp @@ -429,6 +429,27 @@ BOOST_MATH_NAMESPACE_BEGIN return static_cast(-1.5f); } // kurtosis_excess + BOOST_MATH_EXPORT template + BOOST_MATH_GPU_ENABLED inline RealType entropy(const arcsine_distribution& dist) + { + BOOST_MATH_STD_USING + RealType result; + RealType x_min = dist.x_min(); + RealType x_max = dist.x_max(); + + if (!arcsine_detail::check_dist( + "boost::math::entropy(arcsine_distribution<%1%> const&, %1% )", + x_min, + x_max, + &result, Policy()) + ) + { + return result; + } + // A beta(1/2, 1/2) distribution scaled to [x_min, x_max]: log(pi (x_max - x_min) / 4). + return log(x_max - x_min) + log(constants::pi() / 4); + } // entropy + BOOST_MATH_EXPORT template BOOST_MATH_GPU_ENABLED inline RealType kurtosis(const arcsine_distribution& dist) { diff --git a/include/boost/math/distributions/bernoulli.hpp b/include/boost/math/distributions/bernoulli.hpp index f235a62b70..700a6874fe 100644 --- a/include/boost/math/distributions/bernoulli.hpp +++ b/include/boost/math/distributions/bernoulli.hpp @@ -35,6 +35,7 @@ #include // complements #include // error checks #include // isnan. +#include #include #include @@ -335,6 +336,17 @@ BOOST_MATH_NAMESPACE_BEGIN // return (6 * p * p - 6 * p + 1) / (p * (1 - p)) + 3; } + BOOST_MATH_EXPORT template + BOOST_MATH_GPU_ENABLED inline RealType entropy(const bernoulli_distribution& dist) + { + BOOST_MATH_STD_USING + // -p log(p) - (1 - p) log(1 - p), with 0 log(0) = 0: + RealType p = dist.success_fraction(); + if ((p == 0) || (p == 1)) + return 0; + return -p * log(p) - (1 - p) * BOOST_MATH_NAMESPACE::log1p(-p, Policy()); + } + BOOST_MATH_NAMESPACE_END // This include must be at the end, *after* the accessors diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp index 0875cf67a5..a3a85b7512 100644 --- a/include/boost/math/special_functions/gamma.hpp +++ b/include/boost/math/special_functions/gamma.hpp @@ -1874,6 +1874,31 @@ BOOST_MATH_GPU_ENABLED T lgamma_incomplete_imp(T a, T x, const Policy& pol) return log(gamma_q(a, x, pol)); } +// +// log(x^a e^-x / tgamma(a)), the log of regularised_gamma_prefix, for a >= 1. With agh = a + g - 1/2 +// and d = (x - agh) / agh, a log(x) - x - lgamma(a) is exactly +// a log1pmx(d) + x (1/2 - g) / agh + log(agh / e) / 2 - log(lanczos_sum_expG_scaled(a)), +// which has none of the cancellation of the direct form when a and x are large. +// +template +BOOST_MATH_GPU_ENABLED T log_regularised_gamma_prefix(T a, T x, const Policy& pol, const Lanczos&) +{ + BOOST_MATH_STD_USING + if(a < 1) + return a * log(x) - x - BOOST_MATH_NAMESPACE::lgamma(a, pol); + T agh = a + static_cast(Lanczos::g()) - T(0.5); + T d = ((x - a) - static_cast(Lanczos::g()) + T(0.5)) / agh; + return a * BOOST_MATH_NAMESPACE::log1pmx(d, pol) + x * static_cast(0.5 - Lanczos::g()) / agh + + (log(agh) - 1) / 2 - log(Lanczos::lanczos_sum_expG_scaled(a)); +} + +template +BOOST_MATH_GPU_ENABLED T log_regularised_gamma_prefix(T a, T x, const Policy& pol, const lanczos::undefined_lanczos&) +{ + BOOST_MATH_STD_USING + return a * log(x) - x - BOOST_MATH_NAMESPACE::lgamma(a, pol); +} + // Calculate log of incomplete gamma function template BOOST_MATH_GPU_ENABLED T lgamma_incomplete_lower_imp(T a, T x, const Policy& pol) @@ -1910,7 +1935,21 @@ BOOST_MATH_GPU_ENABLED T lgamma_incomplete_lower_imp(T a, T x, const Policy& pol need_p = true; if (need_p) - return log(gamma_p(a, x, pol)); + { + // Method 2 in logs, P = prefix * series / a. It works where P underflows (#1564), and for + // large a it is more accurate than log(P), since log P is far better conditioned than P + // there. For smaller a, or without a Lanczos approximation for the prefix, log(P) is the + // more accurate of the two. + typedef typename lanczos::lanczos::type lanczos_type; + constexpr bool accurate_log_prefix = !BOOST_MATH_NAMESPACE::is_same::value; + if(accurate_log_prefix && (x < a) && (a > 1000)) + return log_regularised_gamma_prefix(a, x, pol, lanczos_type()) + log(detail::lower_gamma_series(a, x, pol) / a); + typedef typename policies::normalise >::type quiet_policy; + T p = gamma_p(a, x, quiet_policy()); + if(p >= tools::min_value()) + return log(p); + return log_regularised_gamma_prefix(a, x, pol, lanczos_type()) + log(detail::lower_gamma_series(a, x, pol) / a); + } return log1p(-gamma_q(a, x, pol), pol); } diff --git a/nanobind/CMakeLists.txt b/nanobind/CMakeLists.txt index 10d57c42fd..a57ee1b34f 100644 --- a/nanobind/CMakeLists.txt +++ b/nanobind/CMakeLists.txt @@ -54,6 +54,8 @@ endif() set(pkg_dir "${CMAKE_CURRENT_BINARY_DIR}/boost/math") set_target_properties(_ext PROPERTIES LIBRARY_OUTPUT_DIRECTORY "${pkg_dir}") -add_custom_command(TARGET _ext POST_BUILD +# The pure-Python files, copied on every build so that edits to them need no relink: +add_custom_target(boost_math_py ALL + COMMAND "${CMAKE_COMMAND}" -E make_directory "${pkg_dir}" COMMAND "${CMAKE_COMMAND}" -E copy_if_different - "${CMAKE_CURRENT_SOURCE_DIR}/boost/math/__init__.py" "${pkg_dir}/__init__.py") + "${CMAKE_CURRENT_SOURCE_DIR}/boost/math/__init__.py" "${CMAKE_CURRENT_SOURCE_DIR}/boost/math/ulps.py" "${pkg_dir}") diff --git a/nanobind/boost/math/ulps.py b/nanobind/boost/math/ulps.py new file mode 100644 index 0000000000..1d4736d1ed --- /dev/null +++ b/nanobind/boost/math/ulps.py @@ -0,0 +1,151 @@ +# Copyright 2026 Nick Thompson +# Distributed under the Boost Software License, Version 1.0. +# https://www.boost.org/LICENSE_1_0.txt +"""ulps plots against an mpmath reference. + + import boost.math as bm, mpmath as mp + from boost.math import ulps + ulps.plot(bm.digamma, mp.digamma, np.linspace(-5, 5, 2000)) + ulps.heatmap(bm.lgamma_p, lambda a, x: mp.log(mp.gammainc(a, 0, x, regularized=True)), + np.logspace(0, 7, 300), np.linspace(0.01, 0.99, 300), y_of=lambda a, r: a * r) + +Each point is classified as one of: + ulps |f - ref| / ulp(ref), both finite; + overflow f overflowed (raised or returned +-inf) where ref is finite: shown cyan; + nan f raised a domain error or returned NaN where ref is finite: shown green; + handled ref isn't a finite double and f agrees (inf, or underflow to 0): shown gray. +""" + +import math + +import mpmath as mp +import numpy as np + +OVERFLOW, NAN, HANDLED = -1.0, -2.0, -3.0 +SPECIAL = {OVERFLOW: ("#2ec4d6", "spurious overflow"), NAN: ("#4caf50", "NaN / domain error"), + HANDLED: ("#8a8a85", "out of range, handled")} +SURFACE, INK, MUTED = "#1a1a19", "#f0efec", "#c3c2b7" + + +def error(f, ref, *args, dps=50): + """The error of f(*args) against the mpmath reference ref(*args): ulps, or a SPECIAL code.""" + with mp.workdps(dps): + try: + r = ref(*[mp.mpf(a) for a in args]) + except Exception: + return math.nan # No reference here; left blank in the plots. + try: + got = f(*args) + except OverflowError: + got = math.inf + except (ValueError, ArithmeticError): + got = math.nan + if isinstance(r, mp.mpc) or not mp.isfinite(r): + return HANDLED if (math.isinf(got) or math.isnan(got)) else float("inf") + rf = float(r) + if math.isinf(rf) or (rf == 0 and r != 0): + # The answer isn't representable: overflow to inf or underflow to 0 is right. + return HANDLED if (got == rf or (rf == 0 and abs(got) <= np.finfo(float).tiny)) else float( + abs(mp.mpf(got) - r) / np.spacing(np.finfo(float).tiny)) + if math.isinf(got): + return OVERFLOW + if math.isnan(got): + return NAN + return float(abs(mp.mpf(got) - r) / np.spacing(abs(rf)) if rf != 0 else abs(got) / np.spacing(0.0)) + + +def _summary(err): + e = err[err >= 0] + missing = int(np.sum(np.isnan(err))) + counts = {name: int(np.sum(err == code)) for code, (_, name) in SPECIAL.items() if np.any(err == code)} + parts = [f"max {e.max():.3g} ulps, median {np.median(e):.2g}" if e.size else "no finite errors"] + parts += [f"{n} {name}" for name, n in counts.items()] + if missing: + parts.append(f"{missing} without a reference") + return "; ".join(parts) + + +def _style(ax, fig): + fig.set_facecolor(SURFACE) + ax.set_facecolor(SURFACE) + ax.tick_params(colors=MUTED) + for s in ax.spines.values(): + s.set_color(MUTED) + + +def _decades(cb): + from matplotlib.ticker import FuncFormatter + cb.ax.yaxis.set_major_formatter(FuncFormatter(lambda v, _: "$10^{%d}$" % round(np.log10(v)))) + cb.ax.tick_params(colors=MUTED) + cb.outline.set_edgecolor(MUTED) + cb.set_label("error in ulps", color=INK) + + +def plot(f, ref, x, title=None, xlabel="x", xscale="linear", out=None, dps=50, reference=None): + """1D: error in ulps of f(x) against ref(x). Returns (fig, err).""" + import matplotlib.pyplot as plt + x = np.asarray(x, dtype=float) + err = np.array([error(f, ref, xi, dps=dps) for xi in x]) + fig, ax = plt.subplots(figsize=(8, 4.5)) + _style(ax, fig) + ok = err >= 0 + ax.scatter(x[ok], np.maximum(err[ok], 0.01), s=3, color="#f5a524", rasterized=True) + top = np.max(err[ok]) if ok.any() else 1 + for code, (color, name) in SPECIAL.items(): + m = err == code + if m.any(): + ax.scatter(x[m], np.full(m.sum(), max(top, 1) * 2), s=6, color=color, label=name, rasterized=True) + ax.set_xscale(xscale) + ax.set_yscale("log") + ax.set_xlabel(xlabel, color=INK) + ax.set_ylabel("error in ulps", color=INK) + ax.set_title(title or getattr(f, "__name__", ""), color=INK) + if any(np.any(err == c) for c in SPECIAL): + ax.legend(frameon=False, labelcolor=INK) + fig.text(0.5, 0.01, f"Reference: {reference or f'mpmath at {dps} digits'}; {_summary(err)}.", ha="center", + color=MUTED, fontsize=8, wrap=True) + fig.tight_layout(rect=(0, 0.03, 1, 1)) + if out: + fig.savefig(out, dpi=300, facecolor=SURFACE) + return fig, err + + +def heatmap(f, ref, x, y, title=None, xlabel="x", ylabel="y", xscale="log", yscale="log", + y_of=None, out=None, dps=50, reference=None): + """2D: error in ulps of f(x, y) over the grid x by y. With y_of, the second argument is + y_of(x, y) instead of y (e.g. y_of=lambda a, r: a * r to plot against x/a). Returns (fig, err).""" + import matplotlib.pyplot as plt + from matplotlib.colors import ListedColormap, LogNorm + x = np.asarray(x, dtype=float) + y = np.asarray(y, dtype=float) + err = np.empty((y.size, x.size)) + for i, yi in enumerate(y): + for j, xj in enumerate(x): + err[i, j] = error(f, ref, xj, y_of(xj, yi) if y_of else yi, dps=dps) + fig, ax = plt.subplots(figsize=(7.2, 6)) + _style(ax, fig) + finite = np.where(err >= 0, np.maximum(err, 0.1), np.nan) + vmax = max(np.nanmax(finite) if np.isfinite(finite).any() else 1, 1) + im = ax.pcolormesh(x, y, finite, norm=LogNorm(vmin=0.1, vmax=vmax), cmap="inferno", shading="nearest", + rasterized=True) + handles = [] + for code, (color, name) in SPECIAL.items(): + m = err == code + if m.any(): + ax.pcolormesh(x, y, np.where(m, 1.0, np.nan), cmap=ListedColormap([color]), shading="nearest", + rasterized=True) + handles.append(plt.Rectangle((0, 0), 1, 1, color=color, label=name)) + ax.set_xscale(xscale) + ax.set_yscale(yscale) + ax.set_xlabel(xlabel, color=INK, fontsize=12) + ax.set_ylabel(ylabel, color=INK, fontsize=12) + ax.set_title(title or getattr(f, "__name__", ""), color=INK, fontsize=12) + if handles: + ax.legend(handles=handles, frameon=False, labelcolor=INK, loc="upper left", fontsize=8) + _decades(fig.colorbar(im, ax=ax)) + fig.text(0.5, 0.01, f"Reference: {reference or f'mpmath at {dps} digits'}; {y.size} x {x.size} grid.\n" + f"{_summary(err)}; errors below 0.1 ulp shown as 0.1.", ha="center", color=MUTED, fontsize=8) + fig.tight_layout(rect=(0, 0.05, 1, 1)) + if out: + fig.savefig(out, dpi=300, facecolor=SURFACE) + return fig, err diff --git a/nanobind/tests/test_ulps.py b/nanobind/tests/test_ulps.py new file mode 100644 index 0000000000..c0761cbec0 --- /dev/null +++ b/nanobind/tests/test_ulps.py @@ -0,0 +1,18 @@ +# Copyright 2026 Nick Thompson +# Distributed under the Boost Software License, Version 1.0. +# https://www.boost.org/LICENSE_1_0.txt +import math + +import mpmath as mp + +import boost.math as bm +from boost.math import ulps + + +def test_error_classification(): + assert ulps.error(bm.tgamma, mp.gamma, 4.5) < 4 + assert ulps.error(lambda x: 7.0, lambda x: mp.mpf(7) + mp.mpf(2) ** -50, 0.0) == 1 # ulp(7) = 2^-50 + assert ulps.error(lambda x: math.inf, lambda x: mp.mpf(1), 0.0) == ulps.OVERFLOW + assert ulps.error(lambda x: math.nan, lambda x: mp.mpf(1), 0.0) == ulps.NAN + assert ulps.error(bm.tgamma, mp.gamma, 200.0) == ulps.HANDLED + assert math.isnan(ulps.error(bm.tgamma, lambda x: 1 / 0, 2.0)) diff --git a/test/test_arcsine.cpp b/test/test_arcsine.cpp index dc2fba919e..5a4b1f11d2 100644 --- a/test/test_arcsine.cpp +++ b/test/test_arcsine.cpp @@ -297,6 +297,10 @@ void test_spots(RealType) BOOST_CHECK_EQUAL(mean(arcsine_01), 0.5); // 1 / (1 + 1) = 1/2 exactly. BOOST_CHECK_EQUAL(median(arcsine_01), 0.5); // 1 / (1 + 1) = 1/2 exactly. BOOST_CHECK_EQUAL(variance(arcsine_01), 0.125); // 1/8 = 0.125 + // entropy = log(pi (x_max - x_min) / 4), as for beta(1/2, 1/2) scaled to [x_min, x_max]: + BOOST_CHECK_CLOSE_FRACTION(entropy(arcsine_01), static_cast(-0.24156447527049044469103689156329442450L), tolerance); + BOOST_CHECK_CLOSE_FRACTION(entropy(arcsine_distribution(-1, 1)), static_cast(0.45158270528945486472619522989488214357L), tolerance); + BOOST_CHECK_CLOSE_FRACTION(entropy(arcsine_distribution(2, 1002)), static_cast(6.6661908037116466073629374724897981983L), tolerance); BOOST_CHECK_CLOSE_FRACTION(standard_deviation(arcsine_01), one_div_root_two() / 2, tolerance); // 1/ sqrt(s) = 0.35355339059327379 BOOST_CHECK_EQUAL(skewness(arcsine_01), 0); // BOOST_CHECK_EQUAL(kurtosis_excess(arcsine_01), -1.5); // -3/2 diff --git a/test/test_bernoulli.cpp b/test/test_bernoulli.cpp index 256a897783..ad5dfb2064 100644 --- a/test/test_bernoulli.cpp +++ b/test/test_bernoulli.cpp @@ -146,6 +146,16 @@ void test_spots(RealType) static_cast(5.11111111111111111111111111111111111111111111L), tolerance); + // entropy = -p log(p) - (1 - p) log(1 - p): + BOOST_CHECK_CLOSE_FRACTION(entropy(bernoulli_distribution(static_cast(0.1L))), + static_cast(0.3250829733914482395065500282238179392356L), tolerance); + BOOST_CHECK_CLOSE_FRACTION(entropy(bernoulli_distribution(static_cast(0.5L))), + static_cast(0.6931471805599453094172321214581765680755L), tolerance); + BOOST_CHECK_CLOSE_FRACTION(entropy(bernoulli_distribution(static_cast(1e-10L))), + static_cast(2.402585092989045684017824788017689207601e-9L), tolerance); + BOOST_CHECK_EQUAL(entropy(bernoulli_distribution(0)), 0); + BOOST_CHECK_EQUAL(entropy(bernoulli_distribution(1)), 0); + #ifndef BOOST_MATH_NO_EXCEPTIONS BOOST_MATH_CHECK_THROW( quantile( diff --git a/test/test_igamma.hpp b/test/test_igamma.hpp index 5c6461bbb2..a43752741c 100644 --- a/test/test_igamma.hpp +++ b/test/test_igamma.hpp @@ -316,6 +316,18 @@ void test_spots(T, const char* name = nullptr) // Check large a, x values. Precision just isn't great here. BOOST_CHECK_CLOSE(::boost::math::lgamma_p(static_cast(1450.25), static_cast(1500.75)), static_cast(-0.09812447528127799786140178403478691390753413399549580160096975713L), tolerance * (std::is_same::value ? 16 : 1)); BOOST_CHECK_CLOSE(::boost::math::lgamma_p(static_cast(2000), static_cast(1900)), static_cast(-4.448523733381445722945397105917814000790587922314824687065050805L), tolerance * gcc_win_mul * (std::is_same::value ? 8 : 1)); + // log P is finite where P underflows (#1564), and for large a it is computed in logs. Without a + // Lanczos approximation the prefix is a log(x) - x - lgamma(a), which cancels ~100-fold here: + { + T log_tolerance = tolerance * (std::is_same >::type, boost::math::lanczos::undefined_lanczos>::value ? 200 : 4); + BOOST_CHECK_CLOSE(::boost::math::lgamma_p(static_cast(1000000), static_cast(500000)), static_cast(-193154.3141086602556152285324604561949543258601759776033022052943L), log_tolerance); + BOOST_CHECK_CLOSE(::boost::math::lgamma_p(static_cast(10000), static_cast(5000)), static_cast(-1936.302975331597581776736325409680160786732605632025420667837257L), log_tolerance); + if (!std::is_same >::type, boost::math::lanczos::undefined_lanczos>::value) + { + BOOST_CHECK_CLOSE(::boost::math::lgamma_p(static_cast(100000), static_cast(90000)), static_cast(-540.4252806833238788932453223532282504424226780607045633973466037L), log_tolerance); + BOOST_CHECK_CLOSE(::boost::math::lgamma_p(static_cast(1000000), static_cast(990000)), static_cast(-53.56704246543242809844092012858294405874181248142883790959132019L), log_tolerance); + } + } // Without a Lanczos approximation, the prefix x^a e^-x / tgamma(a) near x = a was once ~1e4 eps out. Arb: if (std::is_same >::type, boost::math::lanczos::undefined_lanczos>::value) {