From a6a48229aff72e129ae0d4f2cf9443eac41bf09d Mon Sep 17 00:00:00 2001 From: Felix Agene Date: Sun, 26 Jul 2026 09:53:28 -0500 Subject: [PATCH 1/3] fix: stop prec::ulps_eq! degenerating into a 1e-9 absolute comparison `approx`'s `ulps_eq` short-circuits on `abs_diff_eq(epsilon)` before it consults the ULPs distance. The crate's wrapper paired it with `DEFAULT_EPS = 1e-9`, which made the ULPs bound unreachable, so every internal `ulps_eq!(x, y)` was really a 1e-9 absolute comparison. That matters because the macro is used to recognise *exact* parameter values, so anything within 1e-9 of them silently took a degenerate branch: Binomial(p = 1 - 2^-33, n = 100).pmf(99) 0 (true 1.16e-8) Binomial(p = 1 - 2^-33, n = 100).ln_pmf(99) -inf Geometric(p = 1 - 2^-33).max() 1 (true u64::MAX) Geometric(p = 1 - 2^-33).skewness() inf (true 92681.9) Beta(1 + 2^-33, 1 + 2^-33).pdf(0.0) 1 (true 0) digamma(-1 + 2^-33) -inf (true -8.59e9) Adds `DEFAULT_ULPS_EPS = f64::EPSILON` and uses it for `ulps_eq!` and `assert_ulps_eq!`. The exact values still take the degenerate branches - `Binomial(1.0, 5).pmf(5) == 1`, `Geometric(1.0).max() == 1`, `digamma` at the poles is still -inf - all covered by existing tests. Tests use `1 - 2^-33` rather than `1 - 1e-10` so that `1 - p` is exact and the references are not limited by the representation of `p`. Also fixes `Binomial::entropy`, which summed `-p * ln(p)` unguarded and so returned `NaN` once any mass underflowed to zero (`0 * -inf`) - for `p = 0.5` that is every `n` past about 1100. `Categorical::entropy` already filtered zero terms; this brings `Binomial` in line. --- src/distribution/binomial.rs | 39 ++++++++++++++++++++++++++++++++++- src/distribution/geometric.rs | 12 +++++++++++ src/function/gamma.rs | 19 +++++++++++++++++ src/prec.rs | 15 ++++++++++++-- 4 files changed, 82 insertions(+), 3 deletions(-) diff --git a/src/distribution/binomial.rs b/src/distribution/binomial.rs index c9794f0a..315c26ab 100644 --- a/src/distribution/binomial.rs +++ b/src/distribution/binomial.rs @@ -237,7 +237,11 @@ impl Distribution 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); diff --git a/src/distribution/geometric.rs b/src/distribution/geometric.rs index f74f2a79..3a3f9a78 100644 --- a/src/distribution/geometric.rs +++ b/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); diff --git a/src/function/gamma.rs b/src/function/gamma.rs index bff5c45a..771a62d1 100644 --- a/src/function/gamma.rs +++ b/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() { diff --git a/src/prec.rs b/src/prec.rs index ebf1b6d4..5e23b05b 100644 --- a/src/prec.rs +++ b/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; From 97eca50bcfe9dce89e80871dbda283a20a56f329 Mon Sep 17 00:00:00 2001 From: Felix Agene Date: Sun, 26 Jul 2026 09:54:04 -0500 Subject: [PATCH 2/3] fix: correct Geometric::inverse_cdf against its definition `Geometric` overrides `inverse_cdf` with the closed form `ceil(ln1p(-x) / ln1p(-p))`. The quotient is only *approximately* integral at the step boundaries, so `ceil` lands on either side of the true one and the round-trip fails for most `p`: over a sweep of 200 values of `p` and 400 of `k`, 7255 of 16200 `(p, k)` pairs had `inverse_cdf(cdf(k)) != k`. Being an override, this was untouched by the work on the default `inverse_cdf` (#390, #391), which is why #342 stayed open after the explicit implementation it asks for had landed. The closed form is now corrected against the definition being inverted - the least `k` with `cdf(k) >= x`: * ordinary rounding puts it within one step, so the fast path is two or three `cdf` evaluations; * once `x` is within a few ulp of 1, `1 - x` has lost every significant bit and the closed form can be off without bound (`cdf` has a plateau roughly `1/p` wide there - at `p = 1e-9` the answer is 36331335444), so it falls back to a bisection that is exact by definition. After: 0 genuine round-trip failures (the 402 remaining are all cases where `cdf(k-1) == cdf(k)` in f64, so no integer is recoverable), and 0 violations of the definition over 30266 inputs including `x` within ulps of 0 and 1. This also un-ignores `test_inverse_cdf_small_p` and `test_inverse_cdf_extreme_tail`, marked `#[ignore = "Gaps in pathological corner cases and testing"]`. The former was masking two separate bugs: it needs both this fix and the `ulps_eq!` epsilon fix, because it constructs `p = 1 - 1e-14`, which the old 1e-9 epsilon snapped to the degenerate `p == 1` branch. Cost: 8 -> 23 ns/call for the verifying `cdf` evaluations. Closes #342 --- src/distribution/geometric.rs | 126 +++++++++++++++++++++++++++++----- 1 file changed, 110 insertions(+), 16 deletions(-) diff --git a/src/distribution/geometric.rs b/src/distribution/geometric.rs index 3a3f9a78..2bb4268c 100644 --- a/src/distribution/geometric.rs +++ b/src/distribution/geometric.rs @@ -203,7 +203,48 @@ impl DiscreteCDF 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 } } @@ -437,18 +478,6 @@ 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); @@ -584,6 +613,71 @@ 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); + let mut xs: Vec = (1..400).map(|j| j as f64 / 400.0).collect(); + // x within a few ulp of 1 is where `1 - x` has lost every + // significant bit and the closed form needs the bisection fallback + xs.extend((1..6).map(|u| 1.0 - u as f64 * f64::EPSILON / 2.0)); + xs.extend((1..6).map(|u| u as f64 * f64::MIN_POSITIVE)); + xs.extend((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); @@ -617,8 +711,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)); @@ -631,9 +726,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. From dc1a318dfa9ca19b114c2a6cb5c77f2cdd3da862 Mon Sep 17 00:00:00 2001 From: Felix Agene Date: Sun, 26 Jul 2026 18:45:05 -0500 Subject: [PATCH 3/3] test: build the Geometric inverse_cdf test under no_std The test collected its sample points into `Vec`, which does not exist under `--no-default-features`: the crate is `no_std` without `alloc`. Chaining the four groups of x values gives the same sequence with no allocation. CI did not catch this because the test job runs only default features, and the feature-powerset job passes `--no-dev-deps`, so test code is never compiled without std. Co-Authored-By: Claude Opus 5 (1M context) --- src/distribution/geometric.rs | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/distribution/geometric.rs b/src/distribution/geometric.rs index 2bb4268c..b52e1b19 100644 --- a/src/distribution/geometric.rs +++ b/src/distribution/geometric.rs @@ -644,12 +644,15 @@ mod tests { for i in 1..60 { let p = i as f64 / 60.0; let g = create_ok(p); - let mut xs: Vec = (1..400).map(|j| j as f64 / 400.0).collect(); - // x within a few ulp of 1 is where `1 - x` has lost every - // significant bit and the closed form needs the bisection fallback - xs.extend((1..6).map(|u| 1.0 - u as f64 * f64::EPSILON / 2.0)); - xs.extend((1..6).map(|u| u as f64 * f64::MIN_POSITIVE)); - xs.extend((1..400).map(|j| g.cdf(j))); + // 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 {