From a00e75f9f3dd18e264867f49533e8b11e7c150d8 Mon Sep 17 00:00:00 2001 From: Tony Edgin Date: Fri, 28 Nov 2025 11:07:10 -0700 Subject: [PATCH 01/11] 10889: betaIncomplete NaN handling Bring std.mathspecial.betaIncomplete's handling of NaN parameter values in line with how D's operators like operator!"+" handle them. I.e., the NaN with the largest payload is returned. --- std/internal/math/gammafunction.d | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/std/internal/math/gammafunction.d b/std/internal/math/gammafunction.d index 64c9378bbc8..e80c7cc30b0 100644 --- a/std/internal/math/gammafunction.d +++ b/std/internal/math/gammafunction.d @@ -944,15 +944,21 @@ enum real BETA_BIGINV = 1.084202172485504434007e-19L; */ real betaIncomplete(real aa, real bb, real xx ) { + // If any parameters are NaN, return the NaN with the largest payload. + if (isNaN(aa) || isNaN(bb) || isNaN(xx)) + { + // With cmp, + // -NaN(larger) < -NaN(smaller) < -inf < inf < NaN(smaller) < NaN(larger). + const largerParam = cmp(abs(aa), abs(bb)) >= 0 ? aa : bb; + return cmp(abs(xx), abs(largerParam)) >= 0 ? xx : largerParam; + } + if ( !(aa>0 && bb>0) ) { - if ( isNaN(aa) ) return aa; - if ( isNaN(bb) ) return bb; return real.nan; // domain error } if (!(xx>0 && xx<1.0)) { - if (isNaN(xx)) return xx; if ( xx == 0.0L ) return 0.0; if ( xx == 1.0L ) return 1.0; return real.nan; // domain error From bc524ad3b561a9f88bc805fd84c3895d0fad6986 Mon Sep 17 00:00:00 2001 From: Tony Edgin Date: Fri, 28 Nov 2025 13:30:37 -0700 Subject: [PATCH 02/11] 10889: Added support for degenerate cases to betaIncomplete Support was added to betaIncomplete and betaIncompleteCompl for the cases when a or b is +0 or infinity. --- std/internal/math/gammafunction.d | 39 +++++++++++++++++------- std/mathspecial.d | 50 +++++++++++++++++++++++++++++++ 2 files changed, 79 insertions(+), 10 deletions(-) diff --git a/std/internal/math/gammafunction.d b/std/internal/math/gammafunction.d index e80c7cc30b0..c8a71bc0fef 100644 --- a/std/internal/math/gammafunction.d +++ b/std/internal/math/gammafunction.d @@ -953,16 +953,23 @@ real betaIncomplete(real aa, real bb, real xx ) return cmp(abs(xx), abs(largerParam)) >= 0 ? xx : largerParam; } - if ( !(aa>0 && bb>0) ) - { - return real.nan; // domain error - } - if (!(xx>0 && xx<1.0)) + // domain errors + if (signbit(aa) == 1 || signbit(bb) == 1) return real.nan; + if (xx < 0.0L || xx > 1.0L) return real.nan; + + // edge cases + if ( xx == 0.0L ) return 0.0L; + if ( xx == 1.0L ) return 1.0L; + + // degenerate cases + if (aa is +0.0L || aa is real.infinity || bb is +0.0L || bb is real.infinity) { - if ( xx == 0.0L ) return 0.0; - if ( xx == 1.0L ) return 1.0; - return real.nan; // domain error + if (aa is +0.0L && bb is +0.0L) return real.nan; + if (aa is real.infinity && bb is real.infinity) return real.nan; + if (aa is +0.0L || bb is real.infinity) return 1.0L; + if (aa is real.infinity || bb is +0.0L) return 0.0L; } + if ( (bb * xx) <= 1.0L && xx <= 0.95L) { return betaDistPowerSeries(aa, bb, xx); @@ -1335,11 +1342,23 @@ done: assert(isIdentical(betaIncompleteInv(2,NaN(0xABC),8), NaN(0xABC))); assert(isIdentical(betaIncompleteInv(2,3, NaN(0xABC)), NaN(0xABC))); - assert(isNaN(betaIncomplete(-1, 2, 3))); + assert(isNaN(betaIncomplete(-0., 1, .5))); + assert(isNaN(betaIncomplete(1, -0., .5))); + assert(isNaN(betaIncomplete(1, 1, -1))); + assert(isNaN(betaIncomplete(1, 1, 2))); + + assert(betaIncomplete(+0., +0., 0) == 0); + assert(isNaN(betaIncomplete(+0., +0., .5))); + assert(betaIncomplete(+0., +0., 1) == 1); + assert(betaIncomplete(+0., 1, .5) == 1); + assert(betaIncomplete(1, +0., 0) == 0); + assert(betaIncomplete(1, +0., .5) == 0); + assert(betaIncomplete(1, real.infinity, .5) == 1); + assert(betaIncomplete(real.infinity, real.infinity, 0) == 0); + assert(isNaN(betaIncomplete(real.infinity, real.infinity, .5))); assert(betaIncomplete(1, 2, 0)==0); assert(betaIncomplete(1, 2, 1)==1); - assert(isNaN(betaIncomplete(1, 2, 3))); assert(betaIncompleteInv(1, 1, 0)==0); assert(betaIncompleteInv(1, 1, 1)==1); diff --git a/std/mathspecial.d b/std/mathspecial.d index a26434f6efe..8889cbf6c67 100644 --- a/std/mathspecial.d +++ b/std/mathspecial.d @@ -279,6 +279,19 @@ real logmdigammaInverse(real x) * * Returns: * It returns $(SUB I, x)(a,b), an element of [0,1]. + * + * $(TABLE_SV + * $(TR $(TH a) $(TH b) $(TH x) $(TH betaIncomplete(a, b, x)) ) + * $(TR $(TD negative) $(TD b) $(TD x) $(TD $(NAN)) ) + * $(TR $(TD a) $(TD negative) $(TD x) $(TD $(NAN)) ) + * $(TR $(TD a) $(TD b) $(TD $(LT) 0) $(TD $(NAN)) ) + * $(TR $(TD a) $(TD b) $(TD $(GT) 1) $(TD $(NAN)) ) + * $(TR $(TD +0) $(TD +0) $(TD (0,1)) $(TD $(NAN)) ) + * $(TR $(TD $(INFIN)) $(TD $(INFIN)) $(TD (0,1)) $(TD $(NAN)) ) + * ) + * + * If one or more of the input parameters are $(NAN), the one with the largest payload is returned. + * For equal payloads but with possibly different signs, the order of preference is x, a, b. * * Note: * The integral is evaluated by a continued fraction expansion or, when `b * x` is small, by a @@ -287,10 +300,33 @@ real logmdigammaInverse(real x) * See_Also: $(LREF beta) $(LREF betaIncompleteCompl) */ real betaIncomplete(real a, real b, real x ) +in +{ + if (!isNaN(a) && !isNaN(b) && !isNaN(x)) + { + assert(signbit(a) == 0, "a must be positive"); + assert(signbit(b) == 0, "b must be positive"); + assert(x >= 0 && x <= 1, "x must be in [0,1]"); + } +} +out(i; isNaN(i) || (i >=0 && i <= 1)) +do { return std.internal.math.gammafunction.betaIncomplete(a, b, x); } +/// +@safe unittest +{ + assert(betaIncomplete(1, 1, .5) == .5); + assert(betaIncomplete(+0., +0., 0) == 0); + assert(isNaN(betaIncomplete(+0., +0., .5))); + assert(isNaN(betaIncomplete(real.infinity, real.infinity, .5))); + assert(betaIncomplete(real.infinity, real.infinity, 1) == 1); + assert(betaIncomplete(NaN(0x1), 1, NaN(0x2)) is NaN(0x2)); + assert(betaIncomplete(1, NaN(0x3), -NaN(0x3)) is -NaN(0x3)); +} + /** Regularized incomplete beta function complement $(SUB I$(SUP C), x)(a,b) * * Mathematically, if a $(GT) 0, b $(GT) 0, and 0 $(LE) x $(LE) 1, then @@ -307,6 +343,19 @@ real betaIncomplete(real a, real b, real x ) * * Returns: * It returns $(SUB I$(SUP C), x)(a,b), an element of [0,1]. + * + * $(TABLE_SV + * $(TR $(TH a) $(TH b) $(TH x) $(TH betaIncompleteCompl(a, b, x)) ) + * $(TR $(TD negative) $(TD b) $(TD x) $(TD $(NAN)) ) + * $(TR $(TD a) $(TD negative) $(TD x) $(TD $(NAN)) ) + * $(TR $(TD a) $(TD b) $(TD $(LT) 0) $(TD $(NAN)) ) + * $(TR $(TD a) $(TD b) $(TD $(GT) 1) $(TD $(NAN)) ) + * $(TR $(TD +0) $(TD +0) $(TD (0,1)) $(TD $(NAN)) ) + * $(TR $(TD $(INFIN)) $(TD $(INFIN)) $(TD (0,1)) $(TD $(NAN)) ) + * ) + * + * If one or more of the input parameters are $(NAN), the one with the largest payload is returned. + * For equal payloads but with possibly different signs, the order of preference is x, a, b. * * See_Also: $(LREF beta) $(LREF betaIncomplete) */ @@ -322,6 +371,7 @@ in assert(x >= 0 && x <= 1, "x must be in [0, 1]"); } } +out(i; isNaN(i) || (i >=0 && i <= 1)) do { return std.internal.math.gammafunction.betaIncomplete(b, a, 1-x); From a9c11f86e0b5d393d49bcbc0b4adeae50df95500 Mon Sep 17 00:00:00 2001 From: Tony Edgin Date: Tue, 2 Dec 2025 17:39:34 -0700 Subject: [PATCH 03/11] 10889 betaIncomplete a=b, x=.5 special case --- std/internal/math/gammafunction.d | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/std/internal/math/gammafunction.d b/std/internal/math/gammafunction.d index c8a71bc0fef..d6b4aadbc07 100644 --- a/std/internal/math/gammafunction.d +++ b/std/internal/math/gammafunction.d @@ -970,6 +970,9 @@ real betaIncomplete(real aa, real bb, real xx ) if (aa is real.infinity || bb is +0.0L) return 0.0L; } + // symmetry + if (aa == bb && xx == 0.5L) return 0.5L; + if ( (bb * xx) <= 1.0L && xx <= 0.95L) { return betaDistPowerSeries(aa, bb, xx); @@ -1359,6 +1362,7 @@ done: assert(betaIncomplete(1, 2, 0)==0); assert(betaIncomplete(1, 2, 1)==1); + assert(betaIncomplete(9.99999984824320730e+30, 9.99999984824320730e+30, 0.5) == 0.5L); assert(betaIncompleteInv(1, 1, 0)==0); assert(betaIncompleteInv(1, 1, 1)==1); From 4cda156ca20f323255b3504978905457c99555dd Mon Sep 17 00:00:00 2001 From: Tony Edgin Date: Fri, 5 Dec 2025 17:16:25 -0700 Subject: [PATCH 04/11] 10889: Convert betaIncomplete to use beta --- std/internal/math/gammafunction.d | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/std/internal/math/gammafunction.d b/std/internal/math/gammafunction.d index d6b4aadbc07..38f57527d93 100644 --- a/std/internal/math/gammafunction.d +++ b/std/internal/math/gammafunction.d @@ -1033,7 +1033,7 @@ real betaIncomplete(real aa, real bb, real xx ) t *= pow(x,a); t /= a; t *= w; - t *= gamma(a+b) / (gamma(a) * gamma(b)); + t /= beta(a, b); } else { @@ -1621,8 +1621,7 @@ real betaDistPowerSeries(real a, real b, real x ) u = a * log(x); if ( (a+b) < MAXGAMMA && fabs(u) < MAXLOG ) { - t = gamma(a+b)/(gamma(a)*gamma(b)); - s = s * t * pow(x,a); + s = s * pow(x,a) / beta(a, b); } else { From 39720ba888b37ffe73cf9faa231d91bc22febd0b Mon Sep 17 00:00:00 2001 From: Tony Edgin Date: Sat, 29 Nov 2025 11:18:35 -0700 Subject: [PATCH 05/11] 10889: handle reversal truncation issue in betaIncomplete When x > 1/b and x is larger than the mean, the complement is computed instead. x -> 1-x. If prior to reversal, x < real.epsilon, 1-x = 1, so x = 1 after reversal. Since algorthm requires x < 1, the computation of complement fails. To resolve this failure, when the reversal causes x -> 1, use x -> nextDown(1-x) instead. This fix resolved a problematic case for betaIncompleteInv that is tracked by a unittest. The unittest was updated. --- std/internal/math/gammafunction.d | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/std/internal/math/gammafunction.d b/std/internal/math/gammafunction.d index 38f57527d93..8a91a3018d4 100644 --- a/std/internal/math/gammafunction.d +++ b/std/internal/math/gammafunction.d @@ -992,6 +992,7 @@ real betaIncomplete(real aa, real bb, real xx ) b = aa; xc = xx; x = 1.0L - xx; + if (x == 1.0L) x = nextDown(x); } else { @@ -1363,6 +1364,7 @@ done: assert(betaIncomplete(1, 2, 0)==0); assert(betaIncomplete(1, 2, 1)==1); assert(betaIncomplete(9.99999984824320730e+30, 9.99999984824320730e+30, 0.5) == 0.5L); + assert(betaIncomplete(1.17549435082228751e-38, 9.99999977819630836e+22, 9.99999968265522539e-22) == 1.0L); assert(betaIncompleteInv(1, 1, 0)==0); assert(betaIncompleteInv(1, 1, 1)==1); @@ -1408,8 +1410,11 @@ done: assert(feqrel(betaIncompleteInv(0.01L, b1, 9e-26L), 0x1.549696104490aa9p-830L) >= real.mant_dig-10); // --- Problematic cases --- - // This is a situation where the series expansion fails to converge - assert( isNaN(betaIncompleteInv(0.12167L, 4.0640301659679627772e19L, 0.0813601L))); + // In the past, this was a situation where the series expansion failed + // to converge. + assert(!isNaN(betaIncompleteInv(0.12167L, 4.0640301659679627772e19L, 0.0813601L))); + // Using scipy, the result should be 1.683301919972747e-29. + // This next result is almost certainly erroneous. // Mathematica states: "(cannot be determined by current methods)" assert(betaIncomplete(1.16251e20L, 2.18e39L, 5.45e-20L) == -real.infinity); From e19239a3458d2956f2c8d661ef7bcf35c5a7b7b6 Mon Sep 17 00:00:00 2001 From: Tony Edgin Date: Fri, 5 Dec 2025 18:06:27 -0700 Subject: [PATCH 06/11] 10889: Extracted logBeta into a function --- std/internal/math/gammafunction.d | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/std/internal/math/gammafunction.d b/std/internal/math/gammafunction.d index 8a91a3018d4..9f6a458a6c0 100644 --- a/std/internal/math/gammafunction.d +++ b/std/internal/math/gammafunction.d @@ -878,6 +878,14 @@ real beta(in real x, in real y) } +/* This is the natural logarithm of the absolute value of the beta function. + */ +private real logBeta(real x, in real y) +{ + return logGamma(x) + logGamma(y) - logGamma(x+y); +} + + private { /* * These value can be calculated like this: @@ -1039,7 +1047,7 @@ real betaIncomplete(real aa, real bb, real xx ) else { /* Resort to logarithms. */ - y += t + logGamma(a+b) - logGamma(a) - logGamma(b); + y += t - logBeta(a, b); y += log(w/a); t = exp(y); @@ -1630,7 +1638,8 @@ real betaDistPowerSeries(real a, real b, real x ) } else { - t = logGamma(a+b) - logGamma(a) - logGamma(b) + u + log(s); + t = u + log(s) - logBeta(a, b); + if ( t < MINLOG ) { From a2b848c021c9cfb2d9e5b138af52247f8469646e Mon Sep 17 00:00:00 2001 From: Tony Edgin Date: Fri, 5 Dec 2025 17:28:00 -0700 Subject: [PATCH 07/11] 10889: factor constant out of function for later use --- std/internal/math/gammafunction.d | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/std/internal/math/gammafunction.d b/std/internal/math/gammafunction.d index 9f6a458a6c0..5aa7ee4b5e4 100644 --- a/std/internal/math/gammafunction.d +++ b/std/internal/math/gammafunction.d @@ -422,6 +422,13 @@ real gamma(real x) 4.59084371199880305320475827592915200343410999829340L) >= real.mant_dig-1); } + +/* This is the lower bound on x for when the Stirling approximation can be used + * to compute ln(Γ(x)). + */ +private enum real LN_GAMMA_STIRLING_LB = 13.0L; + + /***************************************************** * Natural logarithm of gamma function. * @@ -480,7 +487,7 @@ real logGamma(real x) return z; } - if ( x < 13.0L ) + if ( x < LN_GAMMA_STIRLING_LB ) { z = 1.0L; nx = floor( x + 0.5L ); From a9c74c454f7c9e40954e7e70a085eb46873394a3 Mon Sep 17 00:00:00 2001 From: Tony Edgin Date: Fri, 5 Dec 2025 18:16:15 -0700 Subject: [PATCH 08/11] 10889: optimized logBeta Combine the Stirling approximations of the individual log gamma terms when one of them is large and the difference between the arguments is significant. This reduces the impact of subtracting the two larger terms. --- std/internal/math/gammafunction.d | 73 ++++++++++++++++++++++++++++--- 1 file changed, 68 insertions(+), 5 deletions(-) diff --git a/std/internal/math/gammafunction.d b/std/internal/math/gammafunction.d index 5aa7ee4b5e4..33aff12f867 100644 --- a/std/internal/math/gammafunction.d +++ b/std/internal/math/gammafunction.d @@ -885,11 +885,64 @@ real beta(in real x, in real y) } -/* This is the natural logarithm of the absolute value of the beta function. +/* This is the natural logarithm of the absolute value of the beta function. It + * tries to eliminate reduce the loss of precision that happens when subtracting + * large numbers by combining the Stirling approximations of the individual + * logGamma calls. + * + * ln|B(x,y)| = ln|Γ(x)| + ln|Γ(y)| - ln|Γ(x+y)|. Stirling's approximation for + * ln|Γ(z)| is ln|Γ(z)| ~ zln(z) - z + ln(2𝜋/z)/2 + 𝚺ₙ₌₁ᴺB₂ₙ/[2n(2n-1)z²ⁿ⁻¹], + * where Bₙ is the nᵗʰ Bernoulli number. + * 𝚺ₙ₌₁ᴺB₂ₙ/[2n(2n-1)z²ⁿ⁻¹] = 𝚺ₙ₌₁ᴺB₂ₙ/[2n(2n-1)z²ⁿ⁻²]/z + * = 𝚺ₙ₌₀ᴺ⁻¹B₂₍ₙ₊₁₎/[(2n+2)(2n+1)z²ⁿ]/z + * = [𝚺ₙ₌₀ᴺ⁻¹Cₙ(1/z²)ⁿ]/z, + * where Cₙ = B₂₍ₙ₊₁₎/[(2n+2)(2n+1)]. + * ln|Γ(z)| ~ zln(z) - z + ln(2𝜋/z)/2 +[𝚺ₙ₌₀ᴺ⁻¹Cₙ(1/z²)ⁿ]/z. */ private real logBeta(real x, in real y) { - return logGamma(x) + logGamma(y) - logGamma(x+y); + const larger = x > y ? x : y; + const smaller = x < y ? x : y; + const sum = larger + smaller; + + if (larger >= LN_GAMMA_STIRLING_LB && sum >= LN_GAMMA_STIRLING_LB && larger - smaller > 10.0L) + { + // Assume x > y + // ln|Γ(x)| - ln|Γ(x+y)| + // ~ x⋅ln(x) - (x+y)ln(x+y) + y + ln(2𝜋/x)/2 - ln(2𝜋/[x+y])/2 + // + [𝚺ₙ₌₀ᴺ⁻¹Cₙ(1/x²)ⁿ]/x - [𝚺ₙ₌₀ᴺ⁻¹Cₙ(1/{x+y}²)ⁿ]/{x+y}. + // x⋅ln(x) - (x+y)ln(x+y) + y + ln(2𝜋/x)/2 - ln(2𝜋/[x+y])/2 + // = ln(xˣ) - ln([x+y]ˣ⁺ʸ) + y + ln(√[2𝜋/x]) - ln(√[2𝜋/{x+y}]) + // = ln(xˣ⁻¹ᐟ²/[x+y]ˣ⁺ʸ⁻¹ᐟ²) + y = ln([x/{x+y}]ˣ⁺ʸ⁻¹ᐟ²x⁻ʸ) + y + // = y - y⋅ln(x) + (.5 - x - y)ln(1 + y/x) + // ln|B(x,y)| + // ~ ln|Γ(y)| + y - y⋅ln(x) + (.5 - x - y)ln(1 + y/x) + // + [𝚺ₙ₌₀ᴺ⁻¹Cₙ(1/x²)ⁿ]/x - [𝚺ₙ₌₀ᴺ⁻¹Cₙ(1/{x+y}²)ⁿ]/{x+y}. + const gamDiffApprox = smaller - smaller*log(larger) + (0.5L - sum)*log1p(smaller/larger); + + const gamDiffCorr + = poly(1.0L/larger^^2, logGammaStirlingCoeffs) / larger + - poly(1.0L/sum^^2, logGammaStirlingCoeffs) / sum; + + return logGamma(smaller) + gamDiffApprox + gamDiffCorr; + } + + return logGamma(smaller) + logGamma(larger) - logGamma(sum); +} + +@safe unittest +{ + assert(isClose(logBeta(1, 1), log(beta(1, 1)))); + assert(isClose(logBeta(3, 2), logBeta(2, 3))); + assert(isClose(exp(logBeta(20, 4)), beta(20, 4))); + assert(isClose(logBeta(30, 40), log(beta(30, 40)))); + + // The following were generated by scipy's betaln function. + assert(feqrel(logBeta(-1.4, -0.4), 1.133_156_234_422_692_6) > double.mant_dig-3); + assert(feqrel(logBeta(-0.5, 1.0), 0.693_147_180_559_945_2) > double.mant_dig-3); + assert(feqrel(logBeta(1.0, 2.0), -0.693_147_180_559_945_3) > double.mant_dig-3); + assert(feqrel(logBeta(14.0, 3.0), -7.426_549_072_397_305) > double.mant_dig-3); + assert(feqrel(logBeta(20.0, 30.0), -33.968_820_791_977_386) > double.mant_dig-3); } @@ -1380,6 +1433,11 @@ done: assert(betaIncomplete(1, 2, 1)==1); assert(betaIncomplete(9.99999984824320730e+30, 9.99999984824320730e+30, 0.5) == 0.5L); assert(betaIncomplete(1.17549435082228751e-38, 9.99999977819630836e+22, 9.99999968265522539e-22) == 1.0L); + assert(betaIncomplete(1.00000001954148138e-25, 1.00000001490116119e-01, 1.17549435082228751e-38) == 1.0L); + assert(isClose( + betaIncomplete(9.99999974737875164e-06, 9.99999998050644787e+18, 9.99999968265522539e-22), + 0.9999596214389047L)); + assert(betaIncompleteInv(1, 1, 0)==0); assert(betaIncompleteInv(1, 1, 1)==1); @@ -1416,9 +1474,14 @@ done: ) > 10); // This next case uncovered a one-bit difference in the FYL2X instruction // between Intel and AMD processors. This difference gets magnified by 2^^38. - // WolframAlpha crashes attempting to calculate this. - assert(feqrel(betaIncompleteInv(0x1.ff1275ae5b939bcap-41L, 4.6713e18L, 0.0813601L), - 0x1.f97749d90c7adba8p-63L) >= real.mant_dig - 39); + // WolframAlpha fails to calculate this. + // scipy says that this is 2.225073858507201e-308 in double precision, + // essentially double.min-normal. + assert(isClose( + betaIncompleteInv(0x1.ff1275ae5b939bcap-41L, 4.6713e18L, 0.0813601L), + 2.225073858507201e-308, + 0, + 1e-40)); real a1 = 3.40483L; assert(betaIncompleteInv(a1, 4.0640301659679627772e19L, 0.545113L) == 0x1.ba8c08108aaf5d14p-109L); real b1 = 2.82847e-25L; From 878d31dcb293c0d44e7fbb0e5ebd77e21adcd533 Mon Sep 17 00:00:00 2001 From: Tony Edgin Date: Fri, 5 Dec 2025 18:25:07 -0700 Subject: [PATCH 09/11] 10889: betaDistPowerSeries handle a*s approx eq 1 In betaDistPowerSeries, when s*a is close to 1, extreme loss of precision can occur when combining terms in log space. This fix resolves this issue. --- std/internal/math/gammafunction.d | 37 +++++++++++++++++++++++++++++-- 1 file changed, 35 insertions(+), 2 deletions(-) diff --git a/std/internal/math/gammafunction.d b/std/internal/math/gammafunction.d index 33aff12f867..873c973628f 100644 --- a/std/internal/math/gammafunction.d +++ b/std/internal/math/gammafunction.d @@ -1434,6 +1434,7 @@ done: assert(betaIncomplete(9.99999984824320730e+30, 9.99999984824320730e+30, 0.5) == 0.5L); assert(betaIncomplete(1.17549435082228751e-38, 9.99999977819630836e+22, 9.99999968265522539e-22) == 1.0L); assert(betaIncomplete(1.00000001954148138e-25, 1.00000001490116119e-01, 1.17549435082228751e-38) == 1.0L); + assert(isClose(betaIncomplete(9.99999983775159024e-18, 9.99999977819630836e+22, 1.00000001954148138e-25), 1.0L)); assert(isClose( betaIncomplete(9.99999974737875164e-06, 9.99999998050644787e+18, 9.99999968265522539e-22), 0.9999596214389047L)); @@ -1482,8 +1483,13 @@ done: 2.225073858507201e-308, 0, 1e-40)); + + // scipy says that this is 8.068764506083944e-20 to double precision. Since this is a + // regression test where the original value isn't a known good value, I' updating the + // test value to the current generated value, which is closer to the scipy value. real a1 = 3.40483L; - assert(betaIncompleteInv(a1, 4.0640301659679627772e19L, 0.545113L) == 0x1.ba8c08108aaf5d14p-109L); + assert(betaIncompleteInv(a1, 4.0640301659679627772e19L, 0.545113L) == 0x1.2a867b1e12b9bdf0p-64L); + real b1 = 2.82847e-25L; assert(feqrel(betaIncompleteInv(0.01L, b1, 9e-26L), 0x1.549696104490aa9p-830L) >= real.mant_dig-10); @@ -1708,8 +1714,35 @@ real betaDistPowerSeries(real a, real b, real x ) } else { - t = u + log(s) - logBeta(a, b); + if (abs(a*s - 1.0L) < 0.01L) + { + // Compute logGamma(a+b) - logGamma(b) + real lnGamma_apb_m_lnGamma_b; + + if (b >= LN_GAMMA_STIRLING_LB) + { + const gamDiffApprox = a - a*log(b) + (0.5L - a - b)*log1p(a/b); + + const gamDiffCorr + = poly(1.0L/b^^2, logGammaStirlingCoeffs) / b + - poly(1.0L/(a+b)^^2, logGammaStirlingCoeffs) / (a+b); + + lnGamma_apb_m_lnGamma_b = -gamDiffApprox - gamDiffCorr; + } + else + { + lnGamma_apb_m_lnGamma_b = logGamma(a+b) - logGamma(b); + } + + // Compute log(s) - logGamma(a) + const ln_s_m_lnGamma_a = log1p(a*s - 1.0L) - log(a) - logGamma(a); + t = lnGamma_apb_m_lnGamma_b + u + ln_s_m_lnGamma_a; + } + else + { + t = u + log(s) - logBeta(a, b); + } if ( t < MINLOG ) { From b974408d1415b1db4c60c35c82cb4af9c55a644f Mon Sep 17 00:00:00 2001 From: Tony Edgin Date: Thu, 4 Dec 2025 17:37:43 -0700 Subject: [PATCH 10/11] 10889: clamp betaDistPowerSeries The Beta CDF has a hard upper limit of 1. Due to round-off error, the result of betaDistPowerSeries can be very slightly larger than 1 when it should be 1. In this case, a value of 1 is returned. To prevent clamping from hiding a more serious bug, if the result would be greater than 1 + 2 epsilon, NaN is returned to indicated the computation failed. --- std/internal/math/gammafunction.d | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/std/internal/math/gammafunction.d b/std/internal/math/gammafunction.d index 873c973628f..967eb23e281 100644 --- a/std/internal/math/gammafunction.d +++ b/std/internal/math/gammafunction.d @@ -1469,9 +1469,11 @@ done: assert(betaIncompleteInv(0.01L, 8e-48L, 9e-26L) == 1-real.epsilon); // Beware: a one-bit change in pow() changes almost all digits in the result! + // scipy says that this is 0.99999_99995_89020_6 (0x1.ffff_fffc_783f_2a7ap-1) + // in double precision. assert(feqrel( betaIncompleteInv(0x1.b3d151fbba0eb18p+1L, 1.2265e-19L, 2.44859e-18L), - 0x1.c0110c8531d0952cp-1L + 0x1.ffff_fffc_783f_2a7ap-1 ) > 10); // This next case uncovered a one-bit difference in the FYL2X instruction // between Intel and AMD processors. This difference gets magnified by 2^^38. @@ -1750,6 +1752,8 @@ real betaDistPowerSeries(real a, real b, real x ) } else s = exp(t); } + + if (s > 1.0L) return (s - 2*real.epsilon <= 1.0L) ? 1.0L : real.nan; return s; } From 8a53f717fadca425d6ce17125494a494c87f0cce Mon Sep 17 00:00:00 2001 From: Tony Edgin Date: Sat, 6 Dec 2025 10:24:54 -0700 Subject: [PATCH 11/11] Rerunning github actions because of timeout