fix: correct Geometric::inverse_cdf against its definition

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

@@ -203,7 +203,48 @@ impl DiscreteCDF<u64, f64> for Geometric {
                 self.p, x
             );
         }
-        k as u64
+        // `ln1p(-x) / ln1p(-p)` is only approximately integral at the step
+        // boundaries, so `ceil` lands on either side of the true one and
+        // `inverse_cdf(cdf(k)) != k` for most `p` (statrs-dev/statrs#342).
+        // Everything below corrects the closed form against the definition
+        // being inverted: the least `k` with `cdf(k) >= x`.
+        let candidate = (k.max(1.0) as u64).max(self.min());
+        let is_answer = |k: u64| self.cdf(k) >= x && (k <= self.min() || self.cdf(k - 1) < x);
+
+        // Ordinary rounding puts the closed form within one step, so this is
+        // the path essentially always taken - two or three `cdf` evaluations.
+        if is_answer(candidate) {
+            return candidate;
+        }
+        if candidate > self.min() && is_answer(candidate - 1) {
+            return candidate - 1;
+        }
+        if is_answer(candidate + 1) {
+            return candidate + 1;
+        }
+
+        // Once `x` is within a few ulp of 1, `1 - x` has lost every significant
+        // bit and the closed form can be off by an unbounded number of steps
+        // (`cdf` has a long plateau there - about 1/p wide). Bisect instead,
+        // which is exact by definition regardless of how far off `candidate` is.
+        let mut hi = candidate;
+        while self.cdf(hi) < x {
+            let Some(doubled) = hi.checked_mul(2) else {
+                hi = u64::MAX;
+                break;
+            };
+            hi = doubled;
+        }
+        let mut lo = self.min();
+        while lo < hi {
+            let mid = lo + (hi - lo) / 2;
+            if self.cdf(mid) >= x {
+                hi = mid;
+            } else {
+                lo = mid + 1;
+            }
+        }
+        lo
     }
 }
 
@@ -572,6 +613,74 @@ mod tests {
         density_util::check_discrete_distribution(&create_ok(1.0), 1);
     }
 
+    /// `inverse_cdf` must return the least `k` with `cdf(k) >= x`, and so must
+    /// round-trip `cdf` wherever `cdf` is injective.
+    ///
+    /// The closed form `ceil(ln1p(-x) / ln1p(-p))` is only approximately
+    /// integral at the step boundaries, so on its own it landed on the wrong
+    /// side for 7255 of 16200 (p, k) pairs (statrs-dev/statrs#342).
+    #[test]
+    fn test_inverse_cdf_round_trips_and_matches_definition() {
+        let mut mismatches = 0;
+        for i in 1..200 {
+            let p = i as f64 / 200.0;
+            let g = create_ok(p);
+            for k in 1..=400u64 {
+                let c = g.cdf(k);
+                if c >= 1.0 {
+                    break;
+                }
+                // only meaningful where `cdf` actually separates k-1 from k
+                if g.cdf(k - 1) < c && g.inverse_cdf(c) != k {
+                    mismatches += 1;
+                }
+            }
+        }
+        assert_eq!(mismatches, 0, "cdf/inverse_cdf round-trip failures");
+    }
+
+    #[test]
+    fn test_inverse_cdf_is_least_k_with_cdf_at_least_x() {
+        for i in 1..60 {
+            let p = i as f64 / 60.0;
+            let g = create_ok(p);
+            // Chained rather than collected, so the test also builds under
+            // `no_std`, where there is no `Vec`. The second and third groups
+            // put x within a few ulp of 1, where `1 - x` has lost every
+            // significant bit and the closed form needs the bisection fallback.
+            let xs = (1..400)
+                .map(|j| j as f64 / 400.0)
+                .chain((1..6).map(|u| 1.0 - u as f64 * f64::EPSILON / 2.0))
+                .chain((1..6).map(|u| u as f64 * f64::MIN_POSITIVE))
+                .chain((1..400).map(|j| g.cdf(j)));
+            for x in xs {
+                // `x == 1` saturates to `u64::MAX` by design, iCDF(1) = +inf
+                if !(0.0..=1.0).contains(&x) || x == 1.0 {
+                    continue;
+                }
+                let k = g.inverse_cdf(x);
+                assert!(g.cdf(k) >= x, "p={p} x={x:e}: cdf({k}) < x");
+                assert!(
+                    k <= g.min() || g.cdf(k - 1) < x,
+                    "p={p} x={x:e}: {k} is not the least such k"
+                );
+            }
+        }
+    }
+
+    /// Tiny `p` gives `cdf` a plateau roughly `1/p` wide near 1, so the
+    /// bisection fallback has to cover an unbounded number of steps.
+    #[test]
+    fn test_inverse_cdf_long_plateau() {
+        for p in [1e-3f64, 1e-6, 1e-9] {
+            let g = create_ok(p);
+            let x = 1.0 - f64::EPSILON / 2.0;
+            let k = g.inverse_cdf(x);
+            assert!(g.cdf(k) >= x, "p={p:e}: cdf({k}) < x");
+            assert!(k <= g.min() || g.cdf(k - 1) < x, "p={p:e}: {k} not least");
+        }
+    }
+
     #[test]
     fn test_inverse_cdf() {
         let invcdf = |arg: f64| move |x: Geometric| x.inverse_cdf(arg);
@@ -605,8 +714,9 @@ mod tests {
         distribution.inverse_cdf(0.5);
     }
 
+    /// Un-ignored: this failed until `inverse_cdf` was corrected against the
+    /// definition (statrs-dev/statrs#342).
     #[test]
-    #[ignore = "Gaps in pathological corner cases and testing"]
     fn test_inverse_cdf_small_p() {
         let ident_prob = |arg: f64| move |x: Geometric| x.cdf(x.inverse_cdf(arg));
         let ident_count = |count: u64| move |x: Geometric| x.inverse_cdf(x.cdf(count));
@@ -619,9 +729,8 @@ mod tests {
         test_relative(0.5, q, ident_prob(q));
     }
 
-    #[test]
-    #[ignore = "Gaps in pathological corner cases and testing"]
     /// large k, small p regime
+    #[test]
     fn test_inverse_cdf_extreme_tail() {
         // k_below(k): probability midpoint of the k-th bucket — strictly between cdf(k-1) and
         // cdf(k) — so inverse_cdf must return k.

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;