fix: stop prec::ulps_eq! degenerating into a 1e-9 absolute comparison

src/distribution/binomial.rs

@@ -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)
@@ -412,6 +416,39 @@ 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());
+    }
+
     #[test]
     fn test_pmf() {
         let pmf = |arg: u64| move |x: Binomial| x.pmf(arg);

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/function/gamma.rs

@@ -1195,6 +1195,25 @@ mod tests {
         assert!(super::checked_gamma_ui(1.0, f64::INFINITY).is_err());
     }
 
+    /// `digamma` has poles at the non-positive integers and returns -inf there.
+    /// The pole test is `prec::ulps_eq!(x.floor(), x)`, whose default epsilon
+    /// was `1e-9` absolute, so a whole neighbourhood around each pole returned
+    /// -inf instead of a large finite value.
+    ///
+    /// The tolerance is loose because `digamma(x < 0)` goes through the
+    /// reflection formula and `(PI * x).tan()` loses roughly
+    /// `ulp(PI) / dist_to_pole` of relative accuracy; what is pinned here is
+    /// that the values are finite and of the right magnitude.
+    #[test]
+    fn test_digamma_near_negative_integer_poles_is_finite() {
+        prec::assert_relative_eq!(digamma(-1.0 + f64::powi(2.0, -33)), -8589934591.5772156646, epsilon = 0.0, max_relative = 1e-5);
+        prec::assert_relative_eq!(digamma(-1.0 - f64::powi(2.0, -33)),  8589934592.4227843348, epsilon = 0.0, max_relative = 1e-5);
+        // the poles themselves are unchanged
+        assert_eq!(digamma(-1.0), f64::NEG_INFINITY);
+        assert_eq!(digamma(0.0), f64::NEG_INFINITY);
+        assert_eq!(digamma(-5.0), f64::NEG_INFINITY);
+    }
+
     // TODO: precision testing could be more accurate
     #[test]
     fn test_digamma() {

src/prec.rs

@@ -29,6 +29,7 @@
 //! - `DEFAULT_RELATIVE_ACC`: 1e-14 for relative comparisons
 //! - `DEFAULT_EPS`: 1e-9 for absolute comparisons
 //! - `DEFAULT_ULPS`: 5 for ULPs comparisons
+//! - `DEFAULT_ULPS_EPS`: `f64::EPSILON`, the absolute floor paired with `DEFAULT_ULPS`
 //!
 //! These defaults should be used unless there is a specific reason to use different
 //! precision levels.
@@ -62,6 +63,16 @@ pub const DEFAULT_EPS: f64 = 1e-9;
 /// Default and target ULPs accuracy for f64 operations
 pub const DEFAULT_ULPS: u32 = 5;
 
+/// Default absolute epsilon for ULPs comparisons.
+///
+/// `approx`'s `ulps_eq` short-circuits on `abs_diff_eq(epsilon)` before it looks
+/// at the ULPs distance, so pairing it with [`DEFAULT_EPS`] (`1e-9`) would make
+/// the ULPs bound unreachable and turn `ulps_eq!(x, y)` into a `1e-9` absolute
+/// comparison. `ulps_eq!` is used inside the crate to recognise exact parameter
+/// values (`p == 1.0`, `x == x.floor()`), so its epsilon has to stay at the
+/// scale of a genuine rounding error.
+pub const DEFAULT_ULPS_EPS: f64 = f64::EPSILON;
+
 /// Compares if two floats are close via `approx::abs_diff_eq`
 /// using a maximum absolute difference (epsilon) of `acc`.
 #[deprecated(since = "0.19.0", note = "Use abs_diff_eq! macro instead")]
@@ -144,7 +155,7 @@ mod macros {
     );
     redefine_two_opt_approx_macro!(
         ulps_eq,
-        { epsilon: crate::prec::DEFAULT_EPS, max_ulps: crate::prec::DEFAULT_ULPS }
+        { epsilon: crate::prec::DEFAULT_ULPS_EPS, max_ulps: crate::prec::DEFAULT_ULPS }
     );
 
     pub(crate) use abs_diff_eq;
@@ -162,7 +173,7 @@ mod macros {
     );
     redefine_two_opt_approx_macro!(
         assert_ulps_eq,
-        { epsilon: crate::prec::DEFAULT_EPS, max_ulps: crate::prec::DEFAULT_ULPS }
+        { epsilon: crate::prec::DEFAULT_ULPS_EPS, max_ulps: crate::prec::DEFAULT_ULPS }
     );
 
     pub(crate) use assert_abs_diff_eq;