From 972015f3d7f8a53656c40aeb7455d16b2e71602e Mon Sep 17 00:00:00 2001 From: Felix Agene Date: Sun, 26 Jul 2026 09:55:37 -0500 Subject: [PATCH 1/3] perf: sample Binomial via BINV/BTPE instead of n Bernoulli trials `Binomial::sample` ran one Bernoulli trial per trial, so a single draw was O(n). The reported case in #183 - 5 categories, n = 100_000, 100 samples via `Multinomial` - took 22.4 ms; numpy was measured at ~400x faster. Replaces it with Kachitvichyanukul & Schmeiser (1988): BINV inversion for `n * min(p, 1-p) < 10`, BTPE rejection from a triangle/parallelogram/two exponential tails envelope otherwise, plus exact shortcuts for p in {0, 1} and a Poisson-limit path when `1 - p` rounds to 1. Binomial(0.4, 1_000).sample() 3.6 us -> 56 ns Binomial(0.4, 100_000).sample() 361 us -> 39 ns Binomial(0.4, 1e9).sample() ~3.6 s -> 39 ns Multinomial(5 cats, n=1e5) x100 22.4 ms -> 22.8 us (981x) `Multinomial` benefits automatically, since it draws conditioned binomials per category - which is what the issue asked for. The implementation is adapted from `rand_distr` (MIT/Apache-2.0, attributed in-source). That matters for one detail: the final acceptance test uses the Stirling-series signs from GSL, which differ from the published paper and were verified by one of the algorithm's original designers. Correctness is checked by a chi-square goodness-of-fit test against the exact pmf over five parameter sets covering every code path (BINV, BTPE, and both flipped variants), with expected-count pooling and a 6-sigma threshold on fixed seeds - so it is deterministic and a failure means a real defect. Plus moment tests for n = 1e9 and the Poisson-limit path, and degenerate-parameter tests. Closes #183 --- src/distribution/binomial.rs | 337 ++++++++++++++++++++++++++++++++++- 1 file changed, 333 insertions(+), 4 deletions(-) diff --git a/src/distribution/binomial.rs b/src/distribution/binomial.rs index c9794f0a..8fe9316b 100644 --- a/src/distribution/binomial.rs +++ b/src/distribution/binomial.rs @@ -115,11 +115,218 @@ impl core::fmt::Display for Binomial { #[cfg(feature = "rand")] #[cfg_attr(docsrs, doc(cfg(feature = "rand")))] impl ::rand::distr::Distribution for Binomial { + /// Samples in `O(1)` expected time (independent of `n`) via the BINV + /// inversion algorithm for `n * min(p, 1 - p) < 10` and the BTPE + /// triangle-parallelogram-exponential rejection algorithm otherwise. + /// + /// Kachitvichyanukul, V. and Schmeiser, B. W. (1988). Binomial random + /// variate generation. Communications of the ACM 31(2), 216-222. + /// fn sample(&self, rng: &mut R) -> u64 { - (0..self.n).fold(0, |acc, _| { - let n: f64 = ::rand::RngExt::random(rng); - if n < self.p { acc + 1 } else { acc } - }) + sample_unchecked(rng, self.n, self.p) + } +} + +// The BINV/BTPE implementation below is adapted from the `rand_distr` crate +// (https://github.com/rust-random/rand_distr, Copyright 2018 Developers of the +// Rand project, Copyright 2016-2017 The Rust Project Developers; MIT or +// Apache-2.0), which implements Kachitvichyanukul & Schmeiser (1988) with the +// corrected Stirling-series signs from GSL. + +/// Samples from a binomial distribution with the given `n` and `p`, without +/// validating the parameters. +#[cfg(feature = "rand")] +pub fn sample_unchecked(rng: &mut R, n: u64, p: f64) -> u64 { + if p <= 0.0 || n == 0 { + return 0; + } + if p >= 1.0 { + return n; + } + + // the distribution is symmetric under p -> 1 - p + let flipped = p > 0.5; + let p = if flipped { 1.0 - p } else { p }; + let q = 1.0 - p; + + if q == 1.0 { + // p is below ~2^-54 (and not flipped, since p <= 0.5): the + // distribution is indistinguishable from Poisson(n * p) to O(p) + return super::poisson::sample_unchecked(rng, n as f64 * p) as u64; + } + + // Threshold for preferring BINV; the paper suggests 10. + let sample = if n as f64 * p < 10.0 { + binv(rng, n, p, q) + } else { + btpe(rng, n, p, q) + }; + if flipped { n - sample } else { sample } +} + +/// BINV: sequential inversion from x = 0. Expected iterations ~ n * p. +#[cfg(feature = "rand")] +fn binv(rng: &mut R, n: u64, p: f64, q: f64) -> u64 { + // BINV can get numerically stuck accumulating `u -= r`; a result beyond + // 110 is > 31 standard deviations out for n * p < 10, so restart instead + // (same guard value as GSL and rand_distr). + const BINV_MAX_X: u64 = 110; + + let s = p / q; + let a = (n as f64 + 1.0) * s; + // q^n, via ln_1p to keep accuracy for small p + let r0 = ((-p).ln_1p() * n as f64).exp(); + + 'restart: loop { + let mut r = r0; + let mut u: f64 = ::rand::RngExt::random(rng); + let mut x = 0u64; + while u > r { + u -= r; + x += 1; + if x > BINV_MAX_X { + continue 'restart; + } + r *= a / (x as f64) - s; + } + return x; + } +} + +/// BTPE: rejection from a triangle + parallelogram + two exponential tails +/// envelope around the scaled pmf. Requires `p <= 0.5` and `n * p >= 10`. +#[cfg(feature = "rand")] +#[allow(clippy::many_single_char_names)] // same names as the reference paper +fn btpe(rng: &mut R, n: u64, p: f64, q: f64) -> u64 { + use core::cmp::Ordering; + + // Below this |y - m| the pmf ratio f(y)/f(m) is evaluated directly. + const SQUEEZE_THRESHOLD: u64 = 20; + + // Step 0: constants as functions of n and p. + let n_f = n as f64; + let np = n_f * p; + let npq = np * q; + let f_m = np + p; + let m = f_m as u64; // mode + // radius (and, height being 1, area) of the triangle region + let p1 = (2.195 * npq.sqrt() - 4.6 * q).floor() + 0.5; + let x_m = m as f64 + 0.5; // tip of the triangle + let x_l = x_m - p1; // left edge of the triangle + let x_r = x_m + p1; // right edge of the triangle + let c = 0.134 + 20.5 / (15.3 + m as f64); + // p1 + area of the parallelogram region + let p2 = p1 * (1.0 + 2.0 * c); + + let lambda = |a: f64| a * (1.0 + 0.5 * a); + let lambda_l = lambda((f_m - x_l) / (f_m - x_l * p)); + let lambda_r = lambda((x_r - f_m) / (x_r * q)); + let p3 = p2 + c / lambda_l; + let p4 = p3 + c / lambda_r; + + loop { + // Step 1: select the region via u; v decides acceptance within it. + let u: f64 = ::rand::RngExt::random::(rng) * p4; + let mut v: f64 = ::rand::RngExt::random(rng); + + let y: u64; + if u <= p1 { + // triangle: accept immediately + return (x_m - p1 * v + u) as u64; + } else if u <= p2 { + // Step 2: parallelogram + let x = x_l + (u - p1) / c; + v = v * c + 1.0 - (x - x_m).abs() / p1; + if v > 1.0 { + continue; + } + y = x as u64; + } else if u <= p3 { + // Step 3: left exponential tail (v == 0 gives -inf and retries) + let y_tmp = x_l + v.ln() / lambda_l; + if y_tmp < 0.0 { + continue; + } + y = y_tmp as u64; + v *= (u - p2) * lambda_l; + } else { + // Step 4: right exponential tail (the `as` cast saturates) + let y_tmp = x_r - v.ln() / lambda_r; + if y_tmp > n_f { + continue; + } + y = y_tmp as u64; + v *= (u - p3) * lambda_r; + } + + // Step 5: acceptance/rejection comparison of v against f(y)/f(m). + let k = y.abs_diff(m); + if k <= SQUEEZE_THRESHOLD || (k as f64) >= 0.5 * npq - 1.0 { + // Step 5.1: evaluate the ratio via the recurrence from the mode. + let s = p / q; + let a = s * (n_f + 1.0); + let mut f = 1.0; + match m.cmp(&y) { + Ordering::Less => { + for i in (m + 1)..=y { + f *= a / (i as f64) - s; + } + } + Ordering::Greater => { + for i in (y + 1)..=m { + f /= a / (i as f64) - s; + } + } + Ordering::Equal => {} + } + if v <= f { + return y; + } + continue; + } + + // Step 5.2: squeeze ln(v) between quadratic bounds on ln(f(y)/f(m)). + let kf = k as f64; + let rho = (kf / npq) * ((kf * (kf / 3.0 + 0.625) + 1.0 / 6.0) / npq + 0.5); + let t = -0.5 * kf * kf / npq; + let alpha = v.ln(); + if alpha < t - rho { + return y; + } + if alpha > t + rho { + continue; + } + + // Step 5.3: final comparison against ln(f(y)/f(m)) via Stirling series. + let x1 = (y + 1) as f64; + let f1 = (m + 1) as f64; + let z = ((n - m) + 1) as f64; + let w = ((n - y) + 1) as f64; + + // 13860/166320 = 1/12, 462/166320 = 1/360, ...: the ln k! tail + let stirling = |a: f64| { + let a2 = a * a; + (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / a2) / a2) / a2) / a2) / a / 166320.0 + }; + + let y_sub_m = if y > m { + (y - m) as f64 + } else { + -((m - y) as f64) + }; + // Sign convention on the Stirling terms follows GSL (verified correct + // by one of the algorithm's original designers), not the paper. + if alpha + <= x_m * (f1 / x1).ln() + + (((n - m) as f64) + 0.5) * (z / w).ln() + + y_sub_m * (w * p / (x1 * q)).ln() + + stirling(f1) + + stirling(z) + - stirling(x1) + - stirling(w) + { + return y; + } } } @@ -567,4 +774,126 @@ mod tests { density_util::check_discrete_distribution(&create_ok(0.3, 5), 5); density_util::check_discrete_distribution(&create_ok(0.7, 10), 10); } + /// Chi-square goodness-of-fit of the BINV/BTPE sampler against the exact + /// pmf, over parameter sets covering every code path: BINV (np < 10), BTPE + /// (np >= 10), and both flipped (p > 0.5) variants. Seeds are fixed, so + /// this is deterministic; the 6-sigma acceptance threshold means a failure + /// indicates a real sampler defect, not chance. + #[cfg(feature = "rand")] + #[test] + fn test_sample_chi_square_goodness_of_fit() { + use ::rand::SeedableRng; + use ::rand::distr::Distribution as _; + use ::rand::rngs::StdRng; + + const SAMPLES: usize = 100_000; + for &(n, p) in &[ + (20u64, 0.3f64), // BINV + (100, 0.4), // BTPE + (1000, 0.02), // BTPE, skewed + (100, 0.93), // BINV, flipped + (2000, 0.995), // BTPE, flipped + ] { + let dist = create_ok(p, n); + let mut rng = StdRng::seed_from_u64(0x5EED + n); + + let mean = n as f64 * p; + let sd = (mean * (1.0 - p)).sqrt(); + let lo = (mean - 6.0 * sd).floor().max(0.0) as u64; + let hi = ((mean + 6.0 * sd).ceil() as u64).min(n); + + let mut counts = vec![0u64; (hi - lo + 1) as usize]; + let mut outside = 0u64; + for _ in 0..SAMPLES { + let x: u64 = dist.sample(&mut rng); + match counts.get_mut((x.max(lo) - lo) as usize) { + Some(c) if x >= lo => *c += 1, + _ => outside += 1, + } + } + // beyond 6 sigma the expected count over 1e5 draws is << 1 + assert!(outside <= 3, "n={n} p={p}: {outside} samples beyond 6 sigma"); + + // pool adjacent cells until each has expected count >= 5 + let mut cells: Vec<(f64, u64)> = Vec::new(); + let mut acc = (0.0_f64, 0u64); + for k in lo..=hi { + acc.0 += SAMPLES as f64 * dist.pmf(k); + acc.1 += counts[(k - lo) as usize]; + if acc.0 >= 5.0 { + cells.push(acc); + acc = (0.0, 0); + } + } + if let Some(last) = cells.last_mut() { + last.0 += acc.0; + last.1 += acc.1; + } + + let chi2: f64 = cells + .iter() + .map(|&(e, o)| (o as f64 - e) * (o as f64 - e) / e) + .sum(); + let df = (cells.len() - 1) as f64; + let threshold = df + 6.0 * (2.0 * df).sqrt(); + assert!( + chi2 < threshold, + "n={n} p={p}: chi2 = {chi2:.1} over {} cells (threshold {threshold:.1})", + cells.len() + ); + } + } + + /// Moment sanity for parameter ranges the chi-square test can't cover: + /// n far too large for the old O(n) Bernoulli sampler, and the + /// Poisson-limit path (q rounds to 1.0, i.e. p below ~2^-54). + #[cfg(feature = "rand")] + #[test] + fn test_sample_extreme_parameters_moments() { + use ::rand::SeedableRng; + use ::rand::distr::Distribution as _; + use ::rand::rngs::StdRng; + + // n = 1e9: a single draw used to cost a billion Bernoulli trials + let dist = create_ok(0.4, 1_000_000_000); + let mut rng = StdRng::seed_from_u64(99); + const SAMPLES: usize = 20_000; + let mean = 4.0e8_f64; + let sd = (1.0e9_f64 * 0.4 * 0.6).sqrt(); + let mut sum = 0.0; + for _ in 0..SAMPLES { + let x: u64 = dist.sample(&mut rng); + assert!((x as f64 - mean).abs() < 8.0 * sd, "sample {x} implausibly far out"); + sum += x as f64; + } + let observed_mean = sum / SAMPLES as f64; + // sample mean of 2e4 draws is within ~6 sd / sqrt(SAMPLES) of the mean + prec::assert_abs_diff_eq!(observed_mean, mean, epsilon = 6.0 * sd / (SAMPLES as f64).sqrt()); + + // Poisson-limit path: p < 2^-54 so that 1 - p rounds to 1.0 + // n * p = 0.01 + let dist = create_ok(1e-18, 10_000_000_000_000_000); + let mut rng = StdRng::seed_from_u64(7); + let mut total = 0u64; + for _ in 0..200_000 { + total += >::sample(&dist, &mut rng); + } + // total ~ Poisson(200_000 * 0.01 = 2000); 6 sigma is +-268 + assert!((total as f64 - 2000.0).abs() < 268.0, "Poisson-limit path total {total}"); + } + + /// Degenerate parameters short-circuit. + #[cfg(feature = "rand")] + #[test] + fn test_sample_degenerate() { + use ::rand::SeedableRng; + use ::rand::rngs::StdRng; + + let mut rng = StdRng::seed_from_u64(3); + for _ in 0..10 { + assert_eq!(>::sample(&create_ok(0.0, 50), &mut rng), 0); + assert_eq!(>::sample(&create_ok(1.0, 50), &mut rng), 50); + assert_eq!(>::sample(&create_ok(0.5, 0), &mut rng), 0); + } + } } From aa03f2dd197e78ccc82f003467879807e1c49146 Mon Sep 17 00:00:00 2001 From: Felix Agene Date: Sun, 26 Jul 2026 18:50:17 -0500 Subject: [PATCH 2/3] test: require std for the BTPE chi-square test The test sizes its histogram and its pooled cells from n and p at run time, so it needs `Vec`, which does not exist under `--no-default-features`: the crate is `no_std` without `alloc`. It was gated on `rand` alone, so `--no-default-features --features rand` failed to compile the test binary. Gated on `std` as well rather than rewritten around fixed-size arrays, since a goodness-of-fit test with dynamic binning genuinely wants allocation. The sampler still has no-std coverage from `test_sample_extreme_parameters_moments`, which uses no collections. 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/binomial.rs | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/distribution/binomial.rs b/src/distribution/binomial.rs index 8fe9316b..a07603b4 100644 --- a/src/distribution/binomial.rs +++ b/src/distribution/binomial.rs @@ -779,7 +779,12 @@ mod tests { /// (np >= 10), and both flipped (p > 0.5) variants. Seeds are fixed, so /// this is deterministic; the 6-sigma acceptance threshold means a failure /// indicates a real sampler defect, not chance. - #[cfg(feature = "rand")] + /// + /// Requires `std` as well as `rand`: the histogram and the pooled cells are + /// sized from `n` and `p` at run time, and the crate is `no_std` without + /// `alloc`. The sampler itself is exercised without `std` by + /// [`test_sample_extreme_parameters_moments`], which needs no collections. + #[cfg(all(feature = "rand", feature = "std"))] #[test] fn test_sample_chi_square_goodness_of_fit() { use ::rand::SeedableRng; From 126777b6758233b9e441ab5dab063e9c7e415aac Mon Sep 17 00:00:00 2001 From: Felix Agene Date: Mon, 27 Jul 2026 03:30:59 -0500 Subject: [PATCH 3/3] fix: address review on the BINV/BTPE sampler Both points from @YeungOnion's review on #409. Reflect the Poisson-limit result. That early return skipped the `flipped` reflection every other path applies. It turns out not to be a live bug -- entering the branch needs `p <= 2^-54`, while reflecting can only produce `p >= 2^-53`, the largest f64 below one being `1 - 2^-53` -- but the margin is a single bit, and the comment that stood there gave the wrong reason for why it was safe. Applying the reflection makes the branch correct on its own terms instead of resting on a spacing argument that a later change to the guard could quietly invalidate. The result is also clamped to `n`, since a Poisson variate has no upper bound where a binomial one does. Two tests come with it: one pinning the one-bit reachability margin so that any future change making the branch reachable fails loudly, and one sampling at `p = 1 - 2^-53` -- the closest reachable point on the flipped side, one bit clear of the guard -- which counts failures rather than successes and so exercises the reflection directly. Build the goodness-of-fit test on `stats_tests::chisquare` rather than the hand-rolled statistic and 6-sigma threshold, which gives a real p-value and reuses the crate's own tested implementation. Binning now covers the whole support instead of a 6-sigma window, which is also what makes the observed and expected totals agree, as `chisquare` requires; the final cell is pinned to the exact remainder because summing the pmf across cells otherwise drifts past its `relative_eq` check. The test was mutation-checked rather than assumed meaningful: biasing BINV's ratio by 0.5% takes it to chi-square 2174 over 15 cells at p = 0. Co-Authored-By: Claude Opus 5 (1M context) --- src/distribution/binomial.rs | 147 ++++++++++++++++++++++++++--------- 1 file changed, 112 insertions(+), 35 deletions(-) diff --git a/src/distribution/binomial.rs b/src/distribution/binomial.rs index a07603b4..a931a676 100644 --- a/src/distribution/binomial.rs +++ b/src/distribution/binomial.rs @@ -150,9 +150,18 @@ pub fn sample_unchecked(rng: &mut R, n: u64, p: f64) -> let q = 1.0 - p; if q == 1.0 { - // p is below ~2^-54 (and not flipped, since p <= 0.5): the - // distribution is indistinguishable from Poisson(n * p) to O(p) - return super::poisson::sample_unchecked(rng, n as f64 * p) as u64; + // `p` is at or below 2^-54, so `1 - p` is not representable and the + // distribution is indistinguishable from Poisson(n * p) to O(p). + // + // Reaching here with `flipped` set would require `p <= 2^-54`, while + // reflection can only produce `p >= 2^-53` (the largest f64 below 1 is + // `1 - 2^-53`), so the reflection below is unreachable. It is applied + // anyway: this branch should be correct on its own terms rather than + // resting on a one-bit spacing argument that a later change to the + // guard could silently invalidate. The clamp is needed for the same + // reason -- a Poisson variate has no upper bound, unlike a binomial one. + let sample = (super::poisson::sample_unchecked(rng, n as f64 * p) as u64).min(n); + return if flipped { n - sample } else { sample }; } // Threshold for preferring BINV; the paper suggests 10. @@ -791,6 +800,8 @@ mod tests { use ::rand::distr::Distribution as _; use ::rand::rngs::StdRng; + use crate::stats_tests::chisquare::chisquare; + const SAMPLES: usize = 100_000; for &(n, p) in &[ (20u64, 0.3f64), // BINV @@ -802,49 +813,62 @@ mod tests { let dist = create_ok(p, n); let mut rng = StdRng::seed_from_u64(0x5EED + n); - let mean = n as f64 * p; - let sd = (mean * (1.0 - p)).sqrt(); - let lo = (mean - 6.0 * sd).floor().max(0.0) as u64; - let hi = ((mean + 6.0 * sd).ceil() as u64).min(n); - - let mut counts = vec![0u64; (hi - lo + 1) as usize]; - let mut outside = 0u64; + let mut counts = vec![0usize; (n + 1) as usize]; for _ in 0..SAMPLES { let x: u64 = dist.sample(&mut rng); - match counts.get_mut((x.max(lo) - lo) as usize) { - Some(c) if x >= lo => *c += 1, - _ => outside += 1, - } + assert!(x <= n, "n={n} p={p}: sampled {x}, outside the support"); + counts[x as usize] += 1; } - // beyond 6 sigma the expected count over 1e5 draws is << 1 - assert!(outside <= 3, "n={n} p={p}: {outside} samples beyond 6 sigma"); - // pool adjacent cells until each has expected count >= 5 - let mut cells: Vec<(f64, u64)> = Vec::new(); - let mut acc = (0.0_f64, 0u64); - for k in lo..=hi { + // Bin over the whole support, pooling adjacent outcomes until each + // cell expects at least 5 -- the usual condition for the chi-square + // approximation. Covering the full support rather than a window + // also makes the observed and expected totals agree exactly, which + // `chisquare` checks. + let mut observed: Vec = Vec::new(); + let mut expected: Vec = Vec::new(); + let mut acc = (0.0f64, 0usize); + for k in 0..=n { acc.0 += SAMPLES as f64 * dist.pmf(k); - acc.1 += counts[(k - lo) as usize]; + acc.1 += counts[k as usize]; if acc.0 >= 5.0 { - cells.push(acc); + expected.push(acc.0); + observed.push(acc.1); acc = (0.0, 0); } } - if let Some(last) = cells.last_mut() { - last.0 += acc.0; - last.1 += acc.1; - } + // Fold the leftover tail into the final cell so nothing is dropped. + *expected.last_mut().unwrap() += acc.0; + *observed.last_mut().unwrap() += acc.1; + + assert_eq!( + observed.iter().sum::(), + SAMPLES, + "n={n} p={p}: binning lost samples" + ); - let chi2: f64 = cells - .iter() - .map(|&(e, o)| (o as f64 - e) * (o as f64 - e) / e) - .sum(); - let df = (cells.len() - 1) as f64; - let threshold = df + 6.0 * (2.0 * df).sqrt(); + // `chisquare` rejects inputs whose observed and expected totals + // disagree, and summing the pmf over the support accumulates enough + // rounding to trip that check. Pin the last cell to the exact + // remainder instead: the sum of the others is unchanged, so the + // total then lands on `SAMPLES` to the last bit (the subtraction is + // exact by Sterbenz, the operands being within a factor of two). + let last = expected.len() - 1; + let rest: f64 = expected[..last].iter().sum(); + expected[last] = SAMPLES as f64 - rest; + debug_assert_eq!(expected.iter().sum::(), SAMPLES as f64); + + let (statistic, pvalue) = chisquare(&observed, Some(&expected), None) + .expect("observed and expected totals agree by construction"); + + // The seeds are fixed, so this is deterministic. The threshold is + // deliberately far out in the tail: a correct sampler yields uniform + // p-values, so 1e-6 will not fire by chance, and anything below it + // indicates a real defect rather than an unlucky run. assert!( - chi2 < threshold, - "n={n} p={p}: chi2 = {chi2:.1} over {} cells (threshold {threshold:.1})", - cells.len() + pvalue > 1e-6, + "n={n} p={p}: chi-square = {statistic:.1} over {} cells, p = {pvalue:.3e}", + observed.len() ); } } @@ -887,6 +911,59 @@ mod tests { assert!((total as f64 - 2000.0).abs() < 268.0, "Poisson-limit path total {total}"); } + /// The Poisson-limit branch reflects its result like every other path, but + /// that reflection is currently unreachable: entering the branch needs + /// `p <= 2^-54`, while reflecting can never produce a `p` below `2^-53`. + /// The margin is a single bit, so pin it -- a later change to either the + /// guard or the reflection that makes the branch reachable should fail here + /// rather than silently return an unreflected variate. + #[test] + fn test_poisson_limit_branch_unreachable_when_flipped() { + let largest_below_one = f64::from_bits(1.0f64.to_bits() - 1); + assert!(largest_below_one > 0.5, "premise: this value gets flipped"); + + // Reflection is exact here (Sterbenz), and bottoms out at 2^-53. + let reflected = 1.0 - largest_below_one; + assert_eq!(reflected, (-53f64).exp2()); + + // The branch guard is `1.0 - p == 1.0`, which needs p <= 2^-54. + assert_eq!(1.0 - (-54f64).exp2(), 1.0, "premise: 2^-54 is the threshold"); + assert_ne!(1.0 - reflected, 1.0, "a flipped p reached the Poisson branch"); + } + + /// Mirror of the Poisson-limit case above, on the flipped side: `p` is the + /// largest value below one, so the reflected `p` is `2^-53` -- one bit clear + /// of the Poisson guard, on the ordinary BINV path. Failures rather than + /// successes are counted, which is what the reflection has to get right. + #[cfg(feature = "rand")] + #[test] + fn test_sample_p_adjacent_to_one() { + use ::rand::SeedableRng; + use ::rand::distr::Distribution as _; + use ::rand::rngs::StdRng; + + let n = 10_000_000_000_000_000u64; + let p = f64::from_bits(1.0f64.to_bits() - 1); + let dist = create_ok(p, n); + let mut rng = StdRng::seed_from_u64(11); + + const DRAWS: usize = 200_000; + let mut failures = 0u64; + for _ in 0..DRAWS { + let x: u64 = dist.sample(&mut rng); + assert!(x <= n, "sample {x} exceeds n"); + failures += n - x; + } + + // failures ~ Poisson(DRAWS * n * 2^-53); 6 sigma either side + let mean = DRAWS as f64 * n as f64 * (-53f64).exp2(); + let tol = 6.0 * mean.sqrt(); + assert!( + (failures as f64 - mean).abs() < tol, + "{failures} failures, expected {mean:.0} +- {tol:.0}" + ); + } + /// Degenerate parameters short-circuit. #[cfg(feature = "rand")] #[test]