Panic in Binomial::sample with small `p` and `n` close to u64::MAX

[mstoeckl]

Summary

I observe that when n is a large fraction of u64::MAX and p is a small floating point value, Binomial::sample sometimes panics with the following error:

thread 'main' panicked at /playground/.cargo/registry/src/index.crates.io-1949cf8c6b5b557f/rand_distr-0.5.1/src/binomial.rs:174:5:
assertion failed: x < (i64::MAX as f64)
stack backtrace:
   0: rust_begin_unwind
             at /rustc/4d91de4e48198da2e33413efdcd9cd2cc0c46688/library/std/src/panicking.rs:692:5
   1: core::panicking::panic_fmt
             at /rustc/4d91de4e48198da2e33413efdcd9cd2cc0c46688/library/core/src/panicking.rs:75:14
   2: core::panicking::panic
             at /rustc/4d91de4e48198da2e33413efdcd9cd2cc0c46688/library/core/src/panicking.rs:145:5
   3: rand_distr::binomial::f64_to_i64
             at ./.cargo/registry/src/index.crates.io-1949cf8c6b5b557f/rand_distr-0.5.1/src/binomial.rs:174:5
   4: rand_distr::binomial::btpe
             at ./.cargo/registry/src/index.crates.io-1949cf8c6b5b557f/rand_distr-0.5.1/src/binomial.rs:344:18
   5: <rand_distr::binomial::Binomial as rand::distr::distribution::Distribution<u64>>::sample
             at ./.cargo/registry/src/index.crates.io-1949cf8c6b5b557f/rand_distr-0.5.1/src/binomial.rs:385:49
   6: playground::main
             at ./src/main.rs:19:17
   7: core::ops::function::FnOnce::call_once
             at ./.rustup/toolchains/stable-x86_64-unknown-linux-gnu/lib/rustlib/src/rust/library/core/src/ops/function.rs:250:5
note: Some details are omitted, run with `RUST_BACKTRACE=full` for a verbose backtrace.

I was using current stable rust (1.85) and rand_distr 0.5.1, rand + rand_chacha v0.9.0, on x86_64 (Ryzen 7 7840 HS). The executable is linked to libm from the Arch Linux build of glibc, version 2.41+r9+ga900dbaf70f0-1.

Code sample:

Starting from a specific random generator state, the issue can be reproduced at the Rust Playground with both Stable and Nightly compilers.

https://play.rust-lang.org/?version=stable&mode=debug&edition=2024&gist=da487cdb127284b69d5e74a2f57c9748

use rand::SeedableRng;
use rand_chacha::ChaCha12Rng;
use rand_distr::{Binomial, Distribution};

fn main() {
    let n: u64 = 13301667360732270876;
    let p: f64 = 6.099468852510934e-7;
    let seed: [u8; 32] = [
        65, 56, 42, 90, 76, 75, 12, 151, 226, 205, 74, 251, 40, 152, 42, 207, 114, 35, 176, 234,
        170, 225, 85, 155, 136, 8, 223, 69, 98, 184, 87, 212,
    ];
    let pos: u128 = 1940220;
    let stream: u64 = 0;

    let mut rng = ChaCha12Rng::from_seed(seed);
    rng.set_stream(stream);
    rng.set_word_pos(pos);
    let bin = Binomial::new(n, p).unwrap();
    let value = bin.sample(&mut rng);
    println!("Output: {}", value);
}

Let me know if you'd like me to provide more information; I will try to debug this further if I ever have enough free time.

Workaround

In case anyone else is affected by this edge case, my current workaround is to use the fact that Binomial(n + m, p) can be sampled by sampling x ~ Binomial(n, p), y ~ Binomial(m, p), and outputting x+y (or adding together more copies of smaller distributions, if necessary; I have yet to experimentally verify an exact upper bound for n below which this error does not occur. n <= 1<<62 appears fine so far.)

[benjamin-lieser]

The issue is in

let z = (f64_to_i64(n) + 1 - m) as f64;
let w = (f64_to_i64(n) - y + 1) as f64;

I do not see a good reason why we would not write

let z = n + 1.0 - m as f64;
let w = n - y as f64 + 1.0;

does someone know why it was written like this?

[dhardy]

@benjamin-lieser you adjusted (just moved I think) those lines in rust-random/rand#1484; they were introduced by @vks in rust-random/rand#794.

As to why the code is like this, I don't know. n was converted from u64 to f64 at the start of the method, and is converted to i64 here. Your (@benjamin-lieser's) suggestions would potentially lose precision but the following may be fine:

let z = n - (m - 1) as f64;
let w = n - (y - 1) as f64;

Note: order of operations is potentially important for floats due to rounding but not for integers unless overflow can occur, which doesn't look likely here.

I don't have the time now, but would be worth checking what those values may be in cases like the one reported (among others).

[mstoeckl]

I've opened an issue (#22) for a related problem, which may complicate any attempts at fixing this. The current binomial distribution sampling code appears to be inaccurate for n>2^53.

[dhardy]

Is there actually cause to bother supporting n>2^53 for the binomial distribution?

[mstoeckl]

I've been using rand_distr::Binomial to implement Vitter's Algorithm D† to efficiently sample and iterate over large (10^6~10^10 element) pseudorandom subsets of integers in {0,..,2^64-1}; this involves sampling from a binomial distribution with random p and n ranging up to u64::MAX. The alternative algorithms do not provide the subset in a streaming fashion, and either require additional memory on top of what is needed to store the set or have comparatively slow (O(n log n) instead of O(n)) post-processing phases.

† J. S. Vitter 1984, "Faster Methods for Random Sampling", https://dl.acm.org/doi/pdf/10.1145/358105.893 . The later paper by D. Ting 2021, "Simple Optimal Algorithms for Random Sampling Without Replacement", https://arxiv.org/abs/2104.05091 mentions that this has a simple implementation if you have a binomial sampler.

[benjamin-lieser]

The issue seems a bit deeper. We use i64 internally which will cause problems with n is larger than this. Currently I am not sure why we cannot use u64. I guess it is anyway a good idea to go over the implementation in detail at some point.