Bug. Infinite loop in Binomial Distribution

[benjamin-lieser]

With the current master branch and the values of
n = 16000000 and
p = 3.1444753148558566e-10
for the Binomial distribution

it got stuck in an infinite loop in the BINV algorithm

I know its not reproducible because you need to get an unlucky random number for it to happen. But it is described in the GSL source code https://github.com/ampl/gsl/blob/47862078af4a204c4220f028656a19ddd3922144/randist/binomial_tpe.c#L100

[vks]

Thanks for the report! I managed to reproduce it with this program:

use rand::SeedableRng;
use rand::distributions::Distribution;

fn main() {
    let dist = rand_distr::Binomial::new(16000000, 3.1444753148558566e-10).unwrap();
    let mut sum: u64 = 0;
    let seed = 2569;
    let mut rng = rand_xoshiro::Xoshiro256PlusPlus::seed_from_u64(seed);
    for _ in 0..100_000 {
        sum = sum.wrapping_add(dist.sample(&mut rng));
    }
    println!("sum: {}", sum);
}

[vks]

Here is a reproduction using the internal test::rng:

use rand::distributions::Distribution;

pub fn rng(seed: u64) -> impl rand::RngCore {
    const INC: u64 = 11634580027462260723;
    rand_pcg::Pcg32::new(seed, INC)
}

fn main() {
    let dist = rand_distr::Binomial::new(16000000, 3.1444753148558566e-10).unwrap();
    let mut sum: u64 = 0;
    let seed = 742;
    let mut rng = rng(seed);
    for _ in 0..100_000 {
        sum = sum.wrapping_add(dist.sample(&mut rng));
    }
    println!("sum: {}", sum);
}