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 <noreply@anthropic.com>

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

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.