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

include/boost/math/distributions/arcsine.hpp

@@ -429,6 +429,27 @@ BOOST_MATH_NAMESPACE_BEGIN
       return static_cast<RealType>(-1.5f);
     } // kurtosis_excess
 
+    BOOST_MATH_EXPORT template <class RealType, class Policy>
+    BOOST_MATH_GPU_ENABLED inline RealType entropy(const arcsine_distribution<RealType, Policy>& 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<RealType>() / 4);
+    } // entropy
+
     BOOST_MATH_EXPORT template <class RealType, class Policy>
     BOOST_MATH_GPU_ENABLED inline RealType kurtosis(const arcsine_distribution<RealType, Policy>& dist)
     {

include/boost/math/distributions/bernoulli.hpp

@@ -35,6 +35,7 @@
 #include <boost/math/distributions/complement.hpp> // complements
 #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks
 #include <boost/math/special_functions/fpclassify.hpp> // isnan.
+#include <boost/math/special_functions/log1p.hpp>
 #include <boost/math/policies/policy.hpp>
 #include <boost/math/policies/error_handling.hpp>
 
@@ -335,6 +336,17 @@ BOOST_MATH_NAMESPACE_BEGIN
       // return (6 * p * p - 6 * p + 1) / (p * (1 - p)) + 3;
     }
 
+    BOOST_MATH_EXPORT template <class RealType, class Policy>
+    BOOST_MATH_GPU_ENABLED inline RealType entropy(const bernoulli_distribution<RealType, Policy>& 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

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 <class T, class Policy, class Lanczos>
+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<T>(Lanczos::g()) - T(0.5);
+   T d = ((x - a) - static_cast<T>(Lanczos::g()) + T(0.5)) / agh;
+   return a * BOOST_MATH_NAMESPACE::log1pmx(d, pol) + x * static_cast<T>(0.5 - Lanczos::g()) / agh
+      + (log(agh) - 1) / 2 - log(Lanczos::lanczos_sum_expG_scaled(a));
+}
+
+template <class T, class Policy>
+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 <class T, class Policy>
 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<T, Policy>::type lanczos_type;
+      constexpr bool accurate_log_prefix = !BOOST_MATH_NAMESPACE::is_same<lanczos_type, lanczos::undefined_lanczos>::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<Policy, policies::underflow_error<policies::ignore_error> >::type quiet_policy;
+      T p = gamma_p(a, x, quiet_policy());
+      if(p >= tools::min_value<T>())
+         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);
 }
 

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}")

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

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))

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<RealType>(-0.24156447527049044469103689156329442450L), tolerance);
+    BOOST_CHECK_CLOSE_FRACTION(entropy(arcsine_distribution<RealType>(-1, 1)), static_cast<RealType>(0.45158270528945486472619522989488214357L), tolerance);
+    BOOST_CHECK_CLOSE_FRACTION(entropy(arcsine_distribution<RealType>(2, 1002)), static_cast<RealType>(6.6661908037116466073629374724897981983L), tolerance);
     BOOST_CHECK_CLOSE_FRACTION(standard_deviation(arcsine_01), one_div_root_two<double>() / 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

test/test_bernoulli.cpp

@@ -146,6 +146,16 @@ void test_spots(RealType)
        static_cast<RealType>(5.11111111111111111111111111111111111111111111L),
        tolerance);
 
+  // entropy = -p log(p) - (1 - p) log(1 - p):
+  BOOST_CHECK_CLOSE_FRACTION(entropy(bernoulli_distribution<RealType>(static_cast<RealType>(0.1L))),
+       static_cast<RealType>(0.3250829733914482395065500282238179392356L), tolerance);
+  BOOST_CHECK_CLOSE_FRACTION(entropy(bernoulli_distribution<RealType>(static_cast<RealType>(0.5L))),
+       static_cast<RealType>(0.6931471805599453094172321214581765680755L), tolerance);
+  BOOST_CHECK_CLOSE_FRACTION(entropy(bernoulli_distribution<RealType>(static_cast<RealType>(1e-10L))),
+       static_cast<RealType>(2.402585092989045684017824788017689207601e-9L), tolerance);
+  BOOST_CHECK_EQUAL(entropy(bernoulli_distribution<RealType>(0)), 0);
+  BOOST_CHECK_EQUAL(entropy(bernoulli_distribution<RealType>(1)), 0);
+
 #ifndef BOOST_MATH_NO_EXCEPTIONS
   BOOST_MATH_CHECK_THROW(
      quantile(

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<T>(1450.25), static_cast<T>(1500.75)), static_cast<T>(-0.09812447528127799786140178403478691390753413399549580160096975713L), tolerance * (std::is_same<T, boost::math::concepts::real_concept>::value ? 16 : 1));
    BOOST_CHECK_CLOSE(::boost::math::lgamma_p(static_cast<T>(2000), static_cast<T>(1900)), static_cast<T>(-4.448523733381445722945397105917814000790587922314824687065050805L), tolerance * gcc_win_mul * (std::is_same<T, boost::math::concepts::real_concept>::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<typename boost::math::lanczos::lanczos<T, boost::math::policies::policy<> >::type, boost::math::lanczos::undefined_lanczos>::value ? 200 : 4);
+      BOOST_CHECK_CLOSE(::boost::math::lgamma_p(static_cast<T>(1000000), static_cast<T>(500000)), static_cast<T>(-193154.3141086602556152285324604561949543258601759776033022052943L), log_tolerance);
+      BOOST_CHECK_CLOSE(::boost::math::lgamma_p(static_cast<T>(10000), static_cast<T>(5000)), static_cast<T>(-1936.302975331597581776736325409680160786732605632025420667837257L), log_tolerance);
+      if (!std::is_same<typename boost::math::lanczos::lanczos<T, boost::math::policies::policy<> >::type, boost::math::lanczos::undefined_lanczos>::value)
+      {
+         BOOST_CHECK_CLOSE(::boost::math::lgamma_p(static_cast<T>(100000), static_cast<T>(90000)), static_cast<T>(-540.4252806833238788932453223532282504424226780607045633973466037L), log_tolerance);
+         BOOST_CHECK_CLOSE(::boost::math::lgamma_p(static_cast<T>(1000000), static_cast<T>(990000)), static_cast<T>(-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<typename boost::math::lanczos::lanczos<T, boost::math::policies::policy<> >::type, boost::math::lanczos::undefined_lanczos>::value)
    {