fix: use Loader's saddle-point form for the gamma/beta prefixes and pmfs

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/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/poisson.rs

@@ -276,7 +276,7 @@ impl Discrete<u64, f64> for Poisson {
     ///
     /// where `λ` is the rate
     fn pmf(&self, x: u64) -> f64 {
-        (-self.lambda + x as f64 * self.lambda.ln() - factorial::ln_factorial(x)).exp()
+        self.ln_pmf(x).exp()
     }
 
     /// Calculates the log probability mass function for the poisson
@@ -291,7 +291,19 @@ impl Discrete<u64, f64> for Poisson {
     ///
     /// where `λ` is the rate
     fn ln_pmf(&self, x: u64) -> f64 {
-        -self.lambda + x as f64 * self.lambda.ln() - factorial::ln_factorial(x)
+        let k = x as f64;
+        if x == 0 {
+            return -self.lambda;
+        }
+        if k < gamma::STIRLING_SERIES_MIN {
+            // `ln_factorial` is exact from its table here and every term is
+            // small, so the direct form loses nothing
+            return -self.lambda + k * self.lambda.ln() - factorial::ln_factorial(x);
+        }
+        // Saddle-point form (Loader 2000): the direct expression is a
+        // difference of terms that each grow like `x ln x` while the result
+        // stays `O(1)`, which cost ~1e-11 relative at lambda = 1e4.
+        -gamma::stirling_delta(k) - gamma::bd0(k, self.lambda) - 0.5 * (f64::consts::TAU * k).ln()
     }
 }
 
@@ -421,6 +433,24 @@ mod tests {
         test_exact(10.8, u64::MAX, max);
     }
 
+    /// Large-lambda pmf: the saddle-point form (`gamma::bd0` /
+    /// `gamma::stirling_delta`) replaced `exp(-lam + x ln lam - ln x!)`, whose
+    /// terms each grow like `x ln x` while the result stays `O(1)`. References
+    /// are mpmath at 40 significant digits; the old form was ~650 ulp median and
+    /// 1.4e7 ulp worst over this range.
+    #[test]
+    fn test_pmf_large_lambda_saddle_point() {
+        let pmf = |arg: u64| move |x: Poisson| x.pmf(arg);
+        test_relative(10000.0, 0.0039893895589628256487, pmf(10000));
+        test_relative(1000000.0, 0.0003989422471562440297, pmf(1000000));
+        test_relative(1000000.0, 0.00024189010120174141723, pmf(1001000));
+        // ln_pmf agrees with ln(pmf) where the pmf is comfortably representable
+        for (lam, k) in [(10000.0f64, 10000u64), (1e6, 1001000)] {
+            let d = create_ok(lam);
+            crate::prec::assert_relative_eq!(d.ln_pmf(k), d.pmf(k).ln(), epsilon = 0.0, max_relative = 1e-14);
+        }
+    }
+
     #[test]
     fn test_pmf() {
         let pmf = |arg: u64| move |x: Poisson| x.pmf(arg);