From d89ca19f8ad9302ef418beb3399353f3a414d6ae Mon Sep 17 00:00:00 2001 From: jmh530 Date: Sun, 18 Oct 2026 21:08:54 -0400 Subject: [PATCH] Preserve small shape precision in gamma tail recovery When x/scale is subnormal, the gamma CDF and CCDF recover the lower probability in logarithmic form. Forming 1 + shape for logGamma loses small shape bits even if logGamma itself is accurate. At shape 1e-8 and the smallest positive double argument, this rounding alone causes about 4.7e-12 relative error in the complementary probability. Extend the existing small-shape approximation to include the quadratic term when shape squared is at most the type's epsilon. The omitted cubic term is negligible relative to shape * log(x/scale) in this region. Keep the change within the subnormal-ratio recovery path. Add independent high-precision CDF and CCDF references for the affected case. AI-Assisted-By: Codex (GPT-6 Astra Medium) --- source/mir/stat/distribution/gamma.d | 23 +++++++++++++++++++---- 1 file changed, 19 insertions(+), 4 deletions(-) diff --git a/source/mir/stat/distribution/gamma.d b/source/mir/stat/distribution/gamma.d index 6dcede9..c979e52 100644 --- a/source/mir/stat/distribution/gamma.d +++ b/source/mir/stat/distribution/gamma.d @@ -906,10 +906,12 @@ T gammaSmallRatioLogCDF(T)(const T x, const T shape, const T scale) pragma(inline, false); import mir.math.common: log; import std.mathspecial: logGamma; - // logGamma(1+a)=-EulerGamma*a+O(a*a). Avoid rounding 1+a to 1 - // when this first-order term matters to a small complementary tail. - const T logNext = shape < T.epsilon - ? -T(0.577215664901532860606512090082402431L) * shape + // Forming 1+a loses small shape bits even with an accurate logGamma. + // For a*a <= epsilon, the quadratic expansion omits O(a^3), at most + // O(epsilon*a), which is negligible relative to a*log(x/scale) here. + const T logNext = shape * shape <= T.epsilon + ? shape * (-T(0.577215664901532860606512090082402431L) + + shape * T(0.822467033424113218236207583323012595L)) : T(logGamma(1 + shape)); return shape * (log(x) - log(scale)) - logNext; } @@ -1039,3 +1041,16 @@ version(mir_stat_test) static foreach (T; AliasSeq!(float, double, real)) assert(gammaPDF(T(1), size_t(1_000_000_000)) == 0); } + +// Preserve small shape bits before evaluating logGamma(1+shape). References +// use the exact binary double inputs and high-precision gamma arithmetic. +version(mir_stat_test) +@safe pure nothrow @nogc unittest +{ + import std.math: nextUp; + import mir.math.common: approxEqual; + assert(approxEqual(gammaCCDF(nextUp(0.0), 1e-8), + 7.438600896118198025322571949534157411e-6, 4e-15, 0.0)); + assert(approxEqual(gammaCDF(nextUp(0.0), 1e-8), + 0.999992561399103881801974677428050466, 4e-15, 0.0)); +}