Preserve precision in beta normalization and logarithmic combinatorics

meson.build

@@ -16,6 +16,7 @@ sources_list = [
     'mir/math/internal/powi',
     'mir/math/internal/log_beta',
     'mir/math/internal/log_binomial',
+    'mir/math/internal/log_gamma',
     'mir/math/internal/xlogy',
     'mir/stat/descriptive/multivariate',
     'mir/stat/descriptive/univariate/package',

source/mir/math/internal/log_beta.d

@@ -1,4 +1,11 @@
 /++
+Internal logarithmic beta normalization for positive finite shapes.
+
+Keeps the direct log-gamma calculation for small normal shapes. Larger shapes
+use a combined Stirling expression, so a small difference between large gamma
+logarithms is not lost before the result is returned. Subnormal shapes use the
+gamma recurrence to avoid an overflowing reciprocal inside logGamma.
+
 License: $(HTTP www.apache.org/licenses/LICENSE-2.0, Apache-2.0)
 
 Authors: John Michael Hall
@@ -10,16 +17,139 @@ Copyright: 2022 Mir Stat Authors.
 module mir.math.internal.log_beta;
 
 import mir.internal.utility: isFloatingPoint;
+import mir.math.internal.log_gamma: logGammaStirlingCorrection;
 
-///
+/++
+Natural logarithm of B(alpha, beta). Intermediate arithmetic retains T's
+precision except for the existing Phobos logGamma return type. Results remain
+subject to rounding and may overflow when the logarithm itself exceeds T's range.
++/
 package(mir)
 @safe pure nothrow @nogc
 T logBeta(T)(const T alpha, const T beta)
     if (isFloatingPoint!T)
 {
     import std.mathspecial: logGamma;
 
-    return logGamma(alpha) + logGamma(beta) - logGamma(alpha + beta);
+    if (alpha < 16 && beta < 16 && alpha >= T.min_normal && beta >= T.min_normal)
+        return logGamma(alpha) + logGamma(beta) - logGamma(alpha + beta);
+    return logBetaFallback(alpha, beta);
+}
+
+private @safe pure nothrow @nogc
+T logBetaFallback(T)(const T alpha, const T beta)
+{
+    // Keep the larger calculation out of the ordinary small-shape call site.
+    pragma(inline, false);
+    import std.math: log, log1p;
+    import std.mathspecial: logGamma;
+
+    // Preserve the previous behavior outside the positive finite domain.
+    if (!(alpha > 0 && beta > 0 && alpha < T.infinity && beta < T.infinity))
+        return logGamma(alpha) + logGamma(beta) - logGamma(alpha + beta);
+
+    const T a = alpha < beta ? alpha : beta;
+    const T b = alpha < beta ? beta : alpha;
+    if (b < 16)
+    {
+        // At least one shape is subnormal. Gamma(z)=Gamma(1+z)/z
+        // moves all gamma arguments away from zero; their logs stay finite.
+        return cast(T) (logGamma(1 + a) + logGamma(1 + b) - logGamma(1 + a + b))
+            + log(a + b) - log(a) - log(b);
+    }
+
+    const T ratio = a / b;
+    const T logRatio = log1p(ratio);
+    const T inverseSum = (1 / b) / (1 + ratio);
+    if (a < 16)
+    {
+        const T smallGamma = a < T.min_normal
+            ? cast(T) logGamma(1 + a) - log(a) : cast(T) logGamma(a);
+        // Combine logGamma(b)-logGamma(a+b) before evaluation. Splitting
+        // (a+b-1/2)*log1p(a/b) avoids forming the possibly overflowing sum.
+        return smallGamma - a * log(b) + (a - b * logRatio)
+            - (a - T(0.5)) * logRatio
+            + logGammaStirlingCorrection(1 / b) - logGammaStirlingCorrection(inverseSum);
+    }
+
+    // Both shapes admit Stirling expansion. a/b cannot underflow here:
+    // a>=16 and b<=T.max. log(a/b) avoids subtracting two large logarithms.
+    return a * log(ratio) - a * logRatio - b * logRatio
+        + T(0.5) * logRatio - T(0.5) * log(a)
+        + T(0.91893853320467274178032973640561764L)
+        + logGammaStirlingCorrection(1 / a) + logGammaStirlingCorrection(1 / b)
+        - logGammaStirlingCorrection(inverseSum);
+}
+
+// Independent references around the Stirling transition and unequal shapes.
+version(mir_stat_test)
+@safe pure nothrow @nogc
+unittest
+{
+    import std.meta: AliasSeq;
+    import mir.math.common: approxEqual;
+    static foreach (T; AliasSeq!(float, double, real))
+    {{
+        // 120-decimal-digit references for exactly representable inputs.
+        const real[3][5] references = [
+            [0.125L, 16L, 1.676289181923070711366857286438451596L],
+            [7.25L, 16L, -14.29950746087314117520127940803930057L],
+            [16L, 16L, -22.29368078563352749923792953137939033L],
+            [32L, 64L, -61.71351891820218055826666411894420236L],
+            [0.5L, 1048576L, -6.359106743465463456323875047957258027L]];
+        foreach (row; references)
+        {
+            const actual = logBeta(T(row[0]), T(row[1]));
+            assert(approxEqual(actual, T(row[2]), 64 * T.epsilon, T(0)));
+            assert(actual == logBeta(T(row[1]), T(row[0])));
+        }
+    }}
+}
+
+// Exact beta identities check huge shapes and subnormal gamma arguments.
+version(mir_stat_test)
+@safe pure nothrow @nogc
+unittest
+{
+    import std.meta: AliasSeq;
+    import std.math: log, nextUp, nextDown;
+    import mir.math.common: approxEqual;
+    static foreach (T; AliasSeq!(float, double, real))
+    {{
+        foreach (b; [nextDown(T(16)), T(16), nextUp(T(16)), T(1e16), T.max / 4])
+        {
+            assert(approxEqual(logBeta(T(1), b), -log(b), 128 * T.epsilon, T(0)));
+            assert(approxEqual(logBeta(T(2), b), -log(b) - log(b + 1),
+                128 * T.epsilon, T(0)));
+        }
+        const T tiny = nextUp(T(0));
+        assert(approxEqual(logBeta(tiny, T(1)), -log(tiny), 64 * T.epsilon, T(0)));
+        // B(t,t) tends to 2/t; the correction is negligible at this scale.
+        assert(approxEqual(logBeta(tiny, tiny), log(T(2)) - log(tiny),
+            64 * T.epsilon, T(0)));
+        assert(approxEqual(logBeta(tiny, T(16)), -log(tiny), 64 * T.epsilon, T(0)));
+        const T large = T.max / 4;
+        assert(approxEqual(logBeta(large, large), -large * (2 * log(T(2))),
+            64 * T.epsilon, T(0)));
+    }}
+}
+
+// Inputs outside the positive finite domain retain the gamma-log fallback.
+version(mir_stat_test)
+@safe pure nothrow @nogc
+unittest
+{
+    import std.meta: AliasSeq;
+    import std.math: isNaN;
+    static foreach (T; AliasSeq!(float, double, real))
+    {
+        assert(logBeta(T(0), T(1)) == T.infinity);
+        assert(logBeta(T(1), T(0)) == T.infinity);
+        assert(isNaN(logBeta(T.nan, T(1))));
+        assert(isNaN(logBeta(T(1), T.nan)));
+        assert(isNaN(logBeta(T.infinity, T(1))));
+        assert(isNaN(logBeta(T(1), T.infinity)));
+    }
 }
 
 version(mir_stat_test)

source/mir/math/internal/log_binomial.d

@@ -14,7 +14,11 @@ import mir.internal.utility: isFloatingPoint;
 
 private enum size_t logFactorialAlternative = 1 << 11;
 
-///
+/++
+Natural logarithm of the product of count consecutive integers starting at start.
+An empty product returns zero. The mathematical endpoint may exceed ulong.max;
+large products are evaluated logarithmically without forming that integer sum.
++/
 T logFactorial(T = double)(ulong count, ulong start = 1)
     if (isFloatingPoint!T)
     in (start, "start must be larger than zero")
@@ -23,16 +27,34 @@ T logFactorial(T = double)(ulong count, ulong start = 1)
     import mir.math.numeric: factorial;
     import std.mathspecial: logGamma;
 
-    if (count + start < logFactorialAlternative) {
+    if (count == 0)
+        return T(0);
+    // Compare without wrapping count+start into the small-product path.
+    if (count < logFactorialAlternative && start < logFactorialAlternative - count) {
         return fp_log!T(factorial(count, start));
-    } else {
-        T output = T(logGamma(count + start)); // normally logGamma(x + 1), but start is included in value so subtract out
-        if (start < logFactorialAlternative) {
-            return output - fp_log!T(factorial(start - 1));
-        } else {
-            return output - T(logGamma(start));
-        }
     }
+    // A small lower endpoint has no large gamma logarithm to cancel.
+    // Phobos already evaluates logGamma in real; convert before adding.
+    if (start < 16)
+        return T(logGamma(cast(real) count + start)) - fp_log!T(factorial(start - 1));
+    return logRisingProduct!T(count, start);
+}
+
+private @safe pure nothrow @nogc
+T logRisingProduct(T)(ulong count, ulong start)
+{
+    pragma(inline, false);
+    import std.math: log, log1p;
+    import mir.math.internal.log_gamma: logGammaStirlingCorrection;
+    const T c = count, s = start;
+    const T ratio = c / s;
+    const T logRatio = log1p(ratio);
+    // Combine logGamma(s+c)-logGamma(s) algebraically. Evaluating the two
+    // gamma logs first can lose the entire answer when c is small beside s.
+    // Both gamma arguments are >=16; reciprocal arithmetic avoids s+c.
+    return c * log(s) + (s * logRatio - c) + (c - T(0.5)) * logRatio
+        + logGammaStirlingCorrection((1 / s) / (1 + ratio))
+        - logGammaStirlingCorrection(1 / s);
 }
 
 ///
@@ -85,6 +107,39 @@ unittest
 
 private enum size_t logBinomialCoefficientAlternative = 1 << 11;
 
+// Preserve short products and avoid unsigned wraparound in large endpoints.
+version(mir_stat_test)
+@safe pure nothrow @nogc
+unittest
+{
+    import std.meta: AliasSeq;
+    import std.math: log;
+    import mir.math.common: approxEqual;
+    static foreach (T; AliasSeq!(float, double, real))
+    {{
+        foreach (start; [1UL, 16UL, 2047UL, 2048UL, 10_000_000_000_000_000UL, ulong.max])
+        {
+            assert(logFactorial!T(0, start) == 0);
+            assert(approxEqual(logFactorial!T(1, start), log(T(start)),
+                128 * T.epsilon, T(0)));
+        }
+        // Independent 120-digit references, with exact integer inputs.
+        assert(approxEqual(logFactorial!T(2048, 16),
+            T(13658.47943226613946080615327315399332L), 64 * T.epsilon, T(0)));
+        assert(approxEqual(logFactorial!T(2048, 2047),
+            T(16405.31054988202626462080640908479586L), 64 * T.epsilon, T(0)));
+        assert(approxEqual(logFactorial!T(2, 10_000_000_000_000_000UL),
+            T(73.68272297580946198857572654989964964L), 64 * T.epsilon, T(0)));
+        assert(approxEqual(logFactorial!T(2, ulong.max),
+            T(88.72283911167299960535150143802232549L), 64 * T.epsilon, T(0)));
+        // These endpoints used to wrap into the iterative factorial path.
+        assert(approxEqual(logFactorial!T(ulong.max),
+            T(799877009219260410589.2105935388033377L), 64 * T.epsilon, T(0)));
+        assert(approxEqual(logFactorial!T(ulong.max, ulong.max),
+            T(825449626509665721883.9553386757654554L), 64 * T.epsilon, T(0)));
+    }}
+}
+
 ///
 T logBinomialCoefficient(T = double)(ulong n, uint k)
     if (isFloatingPoint!T)
@@ -105,7 +160,10 @@ T logBinomialCoefficient(T = double)(ulong n, uint k)
         if (n < logBinomialCoefficientAlternative) {
             return fp_log!T(binomialCoefficient(n, k));
         } else {
-            return logFactorial!T(k, n - k + 1) - logFactorial!T(k);
+            import mir.math.internal.log_beta: logBeta;
+            // C(n,k) = 1 / (k * B(n-k+1,k)). k>=2 makes n-k+1 safe
+            // even at ulong.max; no overflowing n+1 or gamma-log difference.
+            return -logBeta(T(n - k + 1), T(k)) - log(T(k));
         }
     }
 }
@@ -203,3 +261,31 @@ unittest
     assert(logBinomialCoefficient(x, cast(uint) x / 2).approxEqual(logFactorial(x) - logFactorial(x / 2) - logFactorial(x / 2)));
     assert(logBinomialCoefficient(x, cast(uint) x - 250).approxEqual(logBinomialCoefficient(x, 250)));
 }
+
+// Large coefficients use beta normalization without n+1 or factorial subtraction.
+version(mir_stat_test)
+@safe pure nothrow @nogc
+unittest
+{
+    import std.meta: AliasSeq;
+    import std.math: log;
+    import mir.math.common: approxEqual;
+    static foreach (T; AliasSeq!(float, double, real))
+    {{
+        foreach (n; [2047UL, 2048UL, 4096UL, 10_000_000_000_000_000UL, ulong.max])
+        {
+            const T expected = log(T(n)) + log(T(n - 1)) - log(T(2));
+            assert(approxEqual(logBinomialCoefficient!T(n, 2), expected,
+                64 * T.epsilon, T(0)));
+            assert(logBinomialCoefficient!T(n, 0) == 0);
+            // Older Phobos returns real even for float/double log inputs.
+            assert(approxEqual(logBinomialCoefficient!T(n, 1), T(log(T(n))),
+                4 * T.epsilon, T(0)));
+        }
+        assert(logBinomialCoefficient!T(2048, 32) == logBinomialCoefficient!T(2048, 2016));
+        assert(approxEqual(logBinomialCoefficient!T(2048, 32),
+            T(162.1864092813636915658867887771846208L), 64 * T.epsilon, T(0)));
+        assert(approxEqual(logBinomialCoefficient!T(ulong.max, uint.max),
+            T(99560390359.53630308783576904322921299L), 64 * T.epsilon, T(0)));
+    }}
+}

source/mir/math/internal/log_gamma.d

@@ -0,0 +1,25 @@
+/++
+Internal Stirling correction shared by logarithmic beta and rising products.
+
+Callers combine the leading log-gamma terms before evaluation to avoid
+cancellation. This helper supplies only the remaining correction at z >= 16,
+using the reciprocal argument so callers need not form a potentially
+unrepresentable endpoint.
+
+License: $(HTTP www.apache.org/licenses/LICENSE-2.0, Apache-2.0)
+Authors: John Michael Hall
+Copyright: 2026 Mir Stat Authors.
++/
+module mir.math.internal.log_gamma;
+
+// Stirling's correction at z>=16, evaluated with r=1/z. The first omitted
+// term is below 2e-23 at z=16, below double and 80-bit real precision.
+package(mir.math) @safe pure nothrow @nogc
+T logGammaStirlingCorrection(T)(const T r)
+{
+    const T r2 = r * r;
+    return r * (T(1) / 12 + r2 * (-T(1) / 360 + r2 * (T(1) / 1260
+        + r2 * (-T(1) / 1680 + r2 * (T(1) / 1188 + r2 * (-T(691) / 360360
+        + r2 * (T(1) / 156 + r2 * (-T(3617) / 122400
+        + r2 * (T(43867) / 244188)))))))));
+}

source/mir/stat/distribution/beta.d

@@ -618,8 +618,10 @@ unittest
 Computes the beta log probability density function (LPDF).
 
 When both shapes are at least 16, a centered calculation avoids subtracting
-large, nearly equal logarithms near the peak. Smaller shapes retain the
-direct logarithmic formula and its accuracy limits.
+large, nearly equal logarithms near the peak. When only one shape is large,
+the beta normalization is evaluated without subtracting large gamma logarithms.
+Two small normal shapes retain the direct normalization and its accuracy limits;
+subnormal shapes use the gamma recurrence to keep that normalization finite.
 
 Returns the natural logarithm of the density.
 
@@ -671,6 +673,31 @@ unittest
     assert(0.25.betaLPDF(0.5, 4).approxEqual(log(betaPDF(0.25, 0.5, 4))));
 }
 
+// A unit shape exposes normalization loss directly, including near an endpoint.
+version(mir_stat_test)
+@safe pure nothrow @nogc
+unittest
+{
+    import std.meta: AliasSeq;
+    import std.math: log, log1p, nextDown;
+    import mir.math.common: approxEqual;
+    static foreach (T; AliasSeq!(float, double, real))
+    {{
+        foreach (a; [T(16), T(1e10), T(1e16), T.max / 4])
+        {
+            assert(approxEqual(betaLPDF(T(1), a, T(1)), log(a),
+                64 * T.epsilon, T(0)));
+            assert(approxEqual(betaLPDF(T(0), T(1), a), log(a),
+                64 * T.epsilon, T(0)));
+            const T x = nextDown(T(1));
+            // B(a,1)=1/a, so log density is log(a)+(a-1)*log(x).
+            const T expected = log(a) + (a - 1) * log1p(x - 1);
+            assert(approxEqual(betaLPDF(x, a, T(1)), expected,
+                64 * T.epsilon, T(0)));
+        }
+    }}
+}
+
 // Stirling's log-gamma correction, with r=1/z and z>=16. The first omitted
 // term is smaller than 2e-23 at z=16, below double/80-bit real precision.
 private @safe pure nothrow @nogc