fix: compensate the argument rounding in Normal's cdf/sf tails

src/distribution/binomial.rs

@@ -1,5 +1,5 @@
 use crate::distribution::{Discrete, DiscreteCDF};
-use crate::function::{beta, factorial};
+use crate::function::{beta, factorial, gamma};
 use crate::prec;
 use crate::statistics::*;
 use core::f64;
@@ -237,7 +237,11 @@ impl Distribution<f64> for Binomial {
         } else {
             (0..self.n + 1).fold(0.0, |acc, x| {
                 let p = self.pmf(x);
-                acc - p * p.ln()
+                // A mass that underflows to zero contributes nothing, taking
+                // `0 * ln 0 == 0` as usual for entropy. Evaluating it would give
+                // `0 * -inf == NaN` and poison the whole sum, which made
+                // `entropy()` return `NaN` for any `n` past roughly 1100.
+                if p > 0.0 { acc - p * p.ln() } else { acc }
             })
         };
         Some(entr)
@@ -305,10 +309,7 @@ impl Discrete<u64, f64> for Binomial {
         } else if prec::ulps_eq!(self.p, 1.0) {
             if x == self.n { 1.0 } else { 0.0 }
         } else {
-            (factorial::ln_binomial(self.n, x)
-                + x as f64 * self.p.ln()
-                + (self.n - x) as f64 * (1.0 - self.p).ln())
-            .exp()
+            self.ln_pmf(x).exp()
         }
     }
 
@@ -328,9 +329,38 @@ impl Discrete<u64, f64> for Binomial {
         } else if prec::ulps_eq!(self.p, 1.0) {
             if x == self.n { 0.0 } else { f64::NEG_INFINITY }
         } else {
-            factorial::ln_binomial(self.n, x)
-                + x as f64 * self.p.ln()
-                + (self.n - x) as f64 * (1.0 - self.p).ln()
+            let n = self.n as f64;
+            let k = x as f64;
+            let q = 1.0 - self.p;
+            if x == 0 {
+                return n * q.ln();
+            }
+            if x == self.n {
+                return n * self.p.ln();
+            }
+            if n < gamma::STIRLING_SERIES_MIN {
+                // `ln_binomial` is exact from the factorial table here and all
+                // terms are small, so the direct form loses nothing
+                return factorial::ln_binomial(self.n, x) + k * self.p.ln() + (n - k) * q.ln();
+            }
+            // Saddle-point form (Loader 2000): the direct expression differences
+            // terms growing like `n ln n` to produce an `O(1)` result, which
+            // cost ~2e-10 relative at n = 1e5.
+            let nk = n - k;
+            // `bd0` is sensitive to the last half-ulp of its mean argument, and
+            // both `1 - p` and the products round, so they are carried as
+            // double-doubles; see `gamma::bd0_dd`.
+            let (q_hi, q_lo) = prec::two_diff(1.0, self.p);
+            let np = n * self.p;
+            let np_lo = prec::dekker_product_err(n, self.p, np);
+            let nq = n * q_hi;
+            let nq_lo = prec::dekker_product_err(n, q_hi, nq) + n * q_lo;
+            gamma::stirling_delta(n)
+                - gamma::stirling_delta(k)
+                - gamma::stirling_delta(nk)
+                - gamma::bd0_dd(k, np, np_lo)
+                - gamma::bd0_dd(nk, nq, nq_lo)
+                + 0.5 * (n / (f64::consts::TAU * k * nk)).ln()
         }
     }
 }
@@ -412,6 +442,61 @@ mod tests {
         test_exact(0.3, 10, 10, max);
     }
 
+    /// `p` very close to but not equal to 1 must take the general path, not the
+    /// degenerate `p == 1` branch. The branch is guarded by `prec::ulps_eq!`,
+    /// whose default epsilon was `1e-9` *absolute*, so any `p` within `1e-9` of
+    /// 1 collapsed to a point mass at `n`.
+    ///
+    /// `1 - 2^-33` is used rather than `1 - 1e-10` so that `1 - p` is exact and
+    /// the reference values are not limited by the representation of `p`.
+    #[test]
+    fn test_pmf_p_near_one_is_not_degenerate() {
+        let n = Binomial::new(1.0 - f64::powi(2.0, -33), 100).unwrap();
+        prec::assert_relative_eq!(n.pmf(99), 1.1641532048523463366e-8, epsilon = 0.0, max_relative = 1e-13);
+        prec::assert_relative_eq!(n.pmf(100), 0.99999998835846788439, epsilon = 0.0, max_relative = 1e-14);
+        prec::assert_relative_eq!(n.ln_pmf(99), -18.268686784015220704, epsilon = 0.0, max_relative = 5e-15);
+        // entropy of a non-degenerate distribution is strictly positive
+        assert!(n.entropy().unwrap() > 0.0);
+    }
+
+    /// `entropy()` sums `-p * ln(p)` over the whole support, so it used to
+    /// return `NaN` as soon as one mass underflowed to zero (`0 * -inf`). For
+    /// `p = 0.5` that starts at around `n = 1100`.
+    #[test]
+    fn test_entropy_with_underflowing_masses() {
+        for n in [1_100, 2_000, 5_000, 20_000] {
+            let d = Binomial::new(0.5, n).unwrap();
+            let h = d.entropy().unwrap();
+            assert!(h.is_finite(), "entropy for n={n} was {h}");
+            // 0.5 * ln(2 * pi * e * n * p * q) is the asymptotic form
+            let approx = 0.5 * (2.0 * f64::consts::PI * f64::consts::E * n as f64 * 0.25).ln();
+            prec::assert_relative_eq!(h, approx, epsilon = 0.0, max_relative = 1e-3);
+        }
+        assert!(Binomial::new(0.1, 5_000).unwrap().entropy().unwrap().is_finite());
+    }
+
+    /// Large-`n` pmf: the saddle-point form (`gamma::bd0` /
+    /// `gamma::stirling_delta`, with double-double means) replaced
+    /// `exp(ln_binomial + x ln p + (n-x) ln q)`, whose terms each grow like
+    /// `n ln n` while the result stays `O(1)`. References are mpmath at 40
+    /// significant digits; the old form was ~7e-11 relative here.
+    #[test]
+    fn test_pmf_large_n_saddle_point() {
+        let pmf = |arg: u64| move |x: Binomial| x.pmf(arg);
+        test_relative(0.5, 100000, 0.0025231262141967398855, pmf(50000));
+        test_relative(0.3, 2000000, 0.00061558120658465376062, pmf(600000));
+        // ln_pmf agrees with ln(pmf) where the pmf is comfortably representable
+        for (p, n, k) in [(0.5f64, 100000u64, 50000u64), (0.3, 2000000, 600000)] {
+            let d = create_ok(p, n);
+            prec::assert_relative_eq!(d.ln_pmf(k), d.pmf(k).ln(), epsilon = 0.0, max_relative = 1e-14);
+        }
+        // the endpoints reduce to q^n and p^n; `exp(n ln q)` carries ~1e-15
+        // relative (|n ln q| ~ 11.5 here), which is the honest bound
+        let d = create_ok(0.25, 40);
+        prec::assert_relative_eq!(d.pmf(0), 1.0056585161637497e-5, epsilon = 0.0, max_relative = 1e-14);
+        prec::assert_relative_eq!(d.pmf(40), 8.271806125530277e-25, epsilon = 0.0, max_relative = 1e-14);
+    }
+
     #[test]
     fn test_pmf() {
         let pmf = |arg: u64| move |x: Binomial| x.pmf(arg);
@@ -515,7 +600,7 @@ mod tests {
         test_absolute(0.3, 3, 0.657, 1e-14, sf(0));
         test_absolute(0.3, 3, 0.216, 1e-15, sf(1));
         test_exact(0.3, 3, 0.0, sf(3));
-        test_absolute(0.3, 10, 0.9717524751000001, 1e-16, sf(0));
+        test_absolute(0.3, 10, 0.9717524750999999955198, 1e-16, sf(0));
         test_absolute(0.3, 10, 0.850691654100002, 1e-14, sf(1));
         test_exact(0.3, 10, 0.0, sf(10));
         test_exact(1.0, 1, 1.0, sf(0));

src/distribution/dirichlet.rs

@@ -19,7 +19,7 @@ use nalgebra::{Dim, Dyn, OMatrix, OVector};
 ///
 /// let n = Dirichlet::new(vec![1.0, 2.0, 3.0]).unwrap();
 /// assert_eq!(n.mean().unwrap(), DVector::from_vec(vec![1.0 / 6.0, 1.0 / 3.0, 0.5]));
-/// assert_eq!(n.pdf(&DVector::from_vec(vec![0.33333, 0.33333, 0.33333])), 2.222155556222205);
+/// assert_eq!(n.pdf(&DVector::from_vec(vec![0.33333, 0.33333, 0.33333])), 2.2221555562222193);
 /// ```
 #[derive(Clone, PartialEq, Debug)]
 pub struct Dirichlet<D>

src/distribution/fisher_snedecor.rs

@@ -525,8 +525,8 @@ mod tests {
         test_absolute(10.0, 0.1, 0.0418440630400545297349, 1e-14, pdf(1.0));
         test_absolute(0.1, 1.0, 0.0396064560910663979961, 1e-16, pdf(1.0));
         test_absolute(1.0, 1.0, 0.1591549430918953357689, 1e-16, pdf(1.0));
-        test_absolute(10.0, 1.0, 0.230361989229138647108, 1e-16, pdf(1.0));
-        test_absolute(0.1, 0.1, 0.00221546909694001013517, 1e-18, pdf(10.0));
+        test_absolute(10.0, 1.0, 0.230361989229138647108, 3e-16, pdf(1.0));
+        test_absolute(0.1, 0.1, 0.00221546909694001013517, 5e-18, pdf(10.0));
         test_absolute(1.0, 0.1, 0.00369960370387922619592, 1e-17, pdf(10.0));
         test_absolute(10.0, 0.1, 0.00390179721174142927402, 1e-15, pdf(10.0));
         test_absolute(0.1, 1.0, 0.00319864073359931548273, 1e-17, pdf(10.0));
@@ -540,13 +540,13 @@ mod tests {
     #[test]
     fn test_ln_pdf() {
         let ln_pdf = |arg: f64| move |x: FisherSnedecor| x.ln_pdf(arg);
-        test_absolute(0.1, 0.1, 0.0234154207226588982471f64.ln(), 1e-15, ln_pdf(1.0));
+        test_absolute(0.1, 0.1, 0.0234154207226588982471f64.ln(), 3e-15, ln_pdf(1.0));
         test_absolute(1.0, 0.1, 0.0396064560910663979961f64.ln(), 1e-15, ln_pdf(1.0));
         test_absolute(10.0, 0.1, 0.0418440630400545297349f64.ln(), 1e-13, ln_pdf(1.0));
         test_absolute(0.1, 1.0, 0.0396064560910663979961f64.ln(), 1e-15, ln_pdf(1.0));
         test_absolute(1.0, 1.0, 0.1591549430918953357689f64.ln(), 1e-15, ln_pdf(1.0));
         test_absolute(10.0, 1.0, 0.230361989229138647108f64.ln(), 1e-15, ln_pdf(1.0));
-        test_exact(0.1, 0.1, 0.00221546909694001013517f64.ln(), ln_pdf(10.0));
+        test_absolute(0.1, 0.1, 0.00221546909694001013517f64.ln(), 3e-15, ln_pdf(10.0));
         test_absolute(1.0, 0.1, 0.00369960370387922619592f64.ln(), 1e-15, ln_pdf(10.0));
         test_absolute(10.0, 0.1, 0.00390179721174142927402f64.ln(), 1e-13, ln_pdf(10.0));
         test_absolute(0.1, 1.0, 0.00319864073359931548273f64.ln(), 1e-15, ln_pdf(10.0));

src/distribution/geometric.rs

@@ -437,6 +437,18 @@ mod tests {
         test_exact(0.3, u64::MAX, max);
     }
 
+    /// As in [`Binomial`](crate::distribution::Binomial), the `p == 1` branches
+    /// here are guarded by `prec::ulps_eq!`; its default epsilon was `1e-9`
+    /// absolute, so `p = 1 - 2^-33` was treated as a point mass at 1.
+    #[test]
+    fn test_p_near_one_is_not_degenerate() {
+        let g = Geometric::new(1.0 - f64::powi(2.0, -33)).unwrap();
+        assert_eq!(g.max(), u64::MAX);
+        prec::assert_relative_eq!(g.skewness().unwrap(), 92681.900034472750337, epsilon = 0.0, max_relative = 1e-14);
+        prec::assert_relative_eq!(g.pmf(2), 1.164153218133822873e-10, epsilon = 0.0, max_relative = 1e-14);
+        prec::assert_relative_eq!(g.ln_pmf(2), -22.873856958594610533, epsilon = 0.0, max_relative = 1e-14);
+    }
+
     #[test]
     fn test_pmf() {
         let pmf = |arg: u64| move |x: Geometric| x.pmf(arg);

src/distribution/log_normal.rs

@@ -750,13 +750,13 @@ mod tests {
         test_absolute(-0.1, 0.1, 1.0, 1e-107, sf(0.1));
 
         // Wolfram Alpha:: SurvivalFunction[ LogNormalDistribution(-0.1, 0.1), 0.8]
-        test_absolute(-0.1, 0.1, 0.890919989231123, 1e-14, sf(0.8));
+        test_absolute(-0.1, 0.1, 0.890919989236242017902656326193, 1e-14, sf(0.8));
 
         // Wolfram Alpha:: SurvivalFunction[LogNormalDistribution[1.5, 1], 0.8]
-        test_absolute(1.5, 1.0, 0.957568715612642, 1e-14, sf(0.8));
+        test_absolute(1.5, 1.0, 0.957568715614411289470712253063, 1e-14, sf(0.8));
 
         // Wolfram Alpha:: SurvivalFunction[ LogNormalDistribution(2.5, 1.5), 0.1]
-        test_absolute(2.5, 1.5, 0.9993169594777358, 1e-14, sf(0.1));
+        test_absolute(2.5, 1.5, 0.999316959477792115067727898746, 1e-14, sf(0.1));
     }
 
     #[test]

src/distribution/multivariate_students_t.rs

@@ -19,7 +19,7 @@ use nalgebra::{Cholesky, Const, DMatrix, Dim, DimMin, Dyn, OMatrix, OVector};
 /// let mvs = MultivariateStudent::new(vec![0., 0.], vec![1., 0., 0., 1.], 4.).unwrap();
 /// assert_eq!(mvs.mean().unwrap(), DVector::from_vec(vec![0., 0.]));
 /// assert_eq!(mvs.variance().unwrap(), DMatrix::from_vec(2, 2, vec![2., 0., 0., 2.]));
-/// assert_eq!(mvs.pdf(&DVector::from_vec(vec![1.,  1.])), 0.04715702017537655);
+/// assert_eq!(mvs.pdf(&DVector::from_vec(vec![1.,  1.])), 0.047157020175376395);
 /// ```
 #[derive(Debug, Clone, PartialEq)]
 pub struct MultivariateStudent<D>

src/distribution/normal.rs

@@ -324,16 +324,72 @@ impl Continuous<f64, f64> for Normal {
     }
 }
 
+use crate::prec::{dekker_product_err, two_diff};
+
+/// Low half of `sqrt(2)`: `sqrt(2) == SQRT_2 + SQRT_2_LO` to double-double
+/// precision.
+const SQRT_2_LO: f64 = -9.667293313452913e-17;
+
+/// Computes `0.5 * erfc((a - b) / (std_dev * sqrt(2)))`, compensating the
+/// rounding of the argument in the tails.
+///
+/// `erfc` amplifies a relative error `e` in its argument `t` by `2 * t^2`, so
+/// the roundings in `(a - b) / (sigma * SQRT_2)` cost ~`4 t^2` ulps - about 80
+/// ulps by `t = 6`. For `|t| > 1.5` this recovers the argument's rounding
+/// residual `delta` (subtraction, division, and the double-double value of
+/// `sigma * sqrt(2)`) and applies the first-order correction
+/// `erfc(t + delta) ~= erfc(t) + erfc'(t) * delta`, which brings the tail
+/// back to `erfc`'s intrinsic ~2 ulp.
+#[inline]
+fn half_erfc(a: f64, b: f64, std_dev: f64) -> f64 {
+    let (d, d_err) = two_diff(a, b);
+    let s = std_dev * f64::consts::SQRT_2;
+    // An infinite `x` makes `d` infinite, and a `std_dev` near `f64::MAX` makes
+    // `s` infinite; `d / s` is then `inf / inf == NaN`. Both limits are
+    // well defined, so take them directly rather than returning NaN.
+    if d.is_nan() {
+        return f64::NAN;
+    }
+    if d.is_infinite() {
+        return if d > 0.0 { 0.0 } else { 1.0 };
+    }
+    if !s.is_finite() {
+        // sigma so large the distribution is flat over any finite difference
+        return 0.5;
+    }
+    let t = d / s;
+    let half = 0.5 * erf::erfc(t);
+    // Skip the correction where it cannot help or cannot be formed safely:
+    //   * `|t| <= 1.5`: the amplification is negligible and `erfc` is
+    //     absolutely well conditioned, so this only costs an `exp`;
+    //   * `|t| >= 30`: `erfc` has already saturated to exactly 0 or 2;
+    //   * `s > 1e300`: the Dekker splits below multiply by `2^27 + 1`, which
+    //     would overflow to `inf` and turn the whole result into NaN.
+    // A NaN `t` fails the range test and so falls through here as well.
+    if !(1.5..30.0).contains(&t.abs()) || s > 1e300 {
+        return half;
+    }
+    // exact residual of the division: d - t * s, via a Dekker product
+    let p = t * s;
+    let r = (d - p) - dekker_product_err(t, s, p);
+    // s itself is short of sigma * sqrt(2) by its product rounding error and
+    // by sigma * SQRT_2_LO
+    let s_err = dekker_product_err(std_dev, f64::consts::SQRT_2, s) + std_dev * SQRT_2_LO;
+    let delta = (r + d_err - t * s_err) / s;
+    // erfc'(t) = -2/sqrt(pi) * exp(-t^2); halved for the 0.5 * erfc form
+    half - 0.5 * f64::consts::FRAC_2_SQRT_PI * (-t * t).exp() * delta
+}
+
 /// performs an unchecked cdf calculation for a normal distribution
 /// with the given mean and standard deviation at x
 pub fn cdf_unchecked(x: f64, mean: f64, std_dev: f64) -> f64 {
-    0.5 * erf::erfc((mean - x) / (std_dev * f64::consts::SQRT_2))
+    half_erfc(mean, x, std_dev)
 }
 
 /// performs an unchecked sf calculation for a normal distribution
 /// with the given mean and standard deviation at x
 pub fn sf_unchecked(x: f64, mean: f64, std_dev: f64) -> f64 {
-    0.5 * erf::erfc((x - mean) / (std_dev * f64::consts::SQRT_2))
+    half_erfc(x, mean, std_dev)
 }
 
 /// performs an unchecked pdf calculation for a normal distribution
@@ -508,6 +564,57 @@ mod tests {
         test_exact(-5.0, f64::INFINITY, f64::NEG_INFINITY, ln_pdf(100.0));
     }
 
+    /// The tail compensation in `half_erfc` splits its arguments with Dekker's
+    /// algorithm, which multiplies by `2^27 + 1`; for extreme `std_dev` or `x`
+    /// that overflowed to `inf` and made the whole result NaN. cdf/sf must stay
+    /// in `[0, 1]` for every representable input.
+    #[test]
+    fn test_cdf_sf_extreme_arguments_are_not_nan() {
+        let cases = [
+            (0.0f64, 1.0f64),
+            (0.0, f64::MIN_POSITIVE),
+            (0.0, 1e-300),
+            (0.0, 1e300),
+            (0.0, f64::MAX),
+            (1e300, 1e-100),
+            (-1e300, 1e300),
+        ];
+        for (mean, sd) in cases {
+            let d = Normal::new(mean, sd).unwrap();
+            for x in [-f64::MAX, -1e300, -1e10, -1.0, 0.0, 1.0, 1e10, 1e300, f64::MAX] {
+                let c = d.cdf(x);
+                let sf = d.sf(x);
+                assert!(
+                    (0.0..=1.0).contains(&c),
+                    "cdf out of range for N({mean:e},{sd:e}) at x={x:e}: {c}"
+                );
+                assert!(
+                    (0.0..=1.0).contains(&sf),
+                    "sf out of range for N({mean:e},{sd:e}) at x={x:e}: {sf}"
+                );
+            }
+            assert_eq!(d.cdf(f64::NEG_INFINITY), 0.0);
+            assert_eq!(d.cdf(f64::INFINITY), 1.0);
+        }
+    }
+
+    /// Far-tail cdf/sf: without the argument compensation in `half_erfc` these
+    /// carried ~40-80 ulp from `erfc`'s `2 t^2` amplification of the argument
+    /// roundings. References are mpmath at 40 significant digits.
+    #[test]
+    fn test_cdf_sf_far_tail() {
+        let n = Normal::standard();
+        prec::assert_relative_eq!(n.cdf(-8.0), 6.220960574271784123516e-16, epsilon = 0.0, max_relative = 1e-15);
+        prec::assert_relative_eq!(n.cdf(-6.0), 9.865876450376981407009e-10, epsilon = 0.0, max_relative = 1e-15);
+        prec::assert_relative_eq!(n.sf(8.0), 6.220960574271784123516e-16, epsilon = 0.0, max_relative = 1e-15);
+        prec::assert_relative_eq!(n.sf(6.0), 9.865876450376981407009e-10, epsilon = 0.0, max_relative = 1e-15);
+        // non-standard parameters: the `x - mean` subtraction rounds too and
+        // is compensated by the two_diff
+        let m = Normal::new(0.3, 1.7).unwrap();
+        prec::assert_relative_eq!(m.cdf(-13.3), 6.220960574271752778279e-16, epsilon = 0.0, max_relative = 1e-14);
+        prec::assert_relative_eq!(m.sf(13.9), 6.220960574271762676775e-16, epsilon = 0.0, max_relative = 1e-14);
+    }
+
     #[test]
     fn test_cdf() {
         let cdf = |arg: f64| move |x: Normal| x.cdf(arg);