From 12f69f99bf40306008c243a89b44c2f122d30f06 Mon Sep 17 00:00:00 2001 From: vangel76 Date: Sat, 7 Nov 2026 17:20:53 +0100 Subject: [PATCH] BINOM.DIST: keep small-k masses whose complement power underflows The small-k product also needs (1 - p)^(n - k) to stay a normal number; otherwise it returned 0 for BINOM.DIST(15,1100,0.5,FALSE) (exact 2.1367646606291324e-298) and the cumulative form built on it. Those cases now use the logarithmic formula. Reference grid gets k = 10..25 at n = 1100. Co-Authored-By: Claude Sonnet 5.5 --- crates/functions/src/stat/distributions.rs | 9 ++++++--- crates/functions/src/stat/tests.rs | 7 +++++++ 2 files changed, 13 insertions(+), 3 deletions(-) diff --git a/crates/functions/src/stat/distributions.rs b/crates/functions/src/stat/distributions.rs index d2e4e3421..0853675a6 100644 --- a/crates/functions/src/stat/distributions.rs +++ b/crates/functions/src/stat/distributions.rs @@ -201,11 +201,14 @@ fn binom_pmf(k: f64, n: f64, p: f64) -> f64 { return p.powf(n); } // A few successes or a few failures: the product form has no large cancelling terms. - // (Not for astronomically large n or a power of p near underflow: the saddle-point form - // below handles those through logs.) + // (Not for astronomically large n, or when either power underflows: the + // saddle-point form below handles those through logs.) if n < 1e30 && p.powf(k) > 1e-280 { if k < 16.0 { - return choose_small(n, k) * p.powi(k as i32) * pow_one_minus(p, n - k); + let rest = pow_one_minus(p, n - k); + if rest >= f64::MIN_POSITIVE { + return choose_small(n, k) * p.powi(k as i32) * rest; + } } if n - k < 16.0 { return choose_small(n, n - k) * p.powf(k) * pow_one_minus(p, n - k); diff --git a/crates/functions/src/stat/tests.rs b/crates/functions/src/stat/tests.rs index 0eb0de0d1..94b35e1f2 100644 --- a/crates/functions/src/stat/tests.rs +++ b/crates/functions/src/stat/tests.rs @@ -745,6 +745,13 @@ const BINOM_REF: &[(f64, f64, f64, f64, f64)] = &[ (99031.0, 100000.0, 0.99, 0.007885853372194001, 0.8416321140066325), (99094.0, 100000.0, 0.99, 0.00013276363518518057, 0.9988439452168577), (99999.0, 100000.0, 0.99, 0.0, 1.0), + // Small k where (1 - p)^(n - k) underflows although the mass is a normal number. + (10.0, 1100.0, 0.5, 5.050689010687844e-308, 5.097367491037505e-308), + (14.0, 1100.0, 0.5, 2.9513324041838846e-300, 2.989803284040754e-300), + (15.0, 1100.0, 0.5, 2.1367646606291324e-298, 2.16666269346954e-298), + (16.0, 1100.0, 0.5, 1.4489935354891303e-296, 1.4706601624238258e-296), + (20.0, 1100.0, 0.5, 1.711087104904223e-289, 1.7433098849775326e-289), + (25.0, 1100.0, 0.5, 3.906996229178926e-281, 3.999839082741416e-281), ]; // Columns: x, mean, mass, cdf. The last block has x > 1000, which takes the incomplete-gamma branch.