Fix u64 overflow in Hypergeometric::new for populations near u64::MAX

CHANGELOG.md

@@ -10,6 +10,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Added a parameter fuzzing script to search for panics and invalid output ([#53])
 
 ### Fixes
+- Fix overflow in `Hypergeometric::new` for populations near `u64::MAX` that corrupted the sampler and could panic ([#60])
 - Avoid returning negative values in `InverseGaussian::sample`; this is a Value-breaking change ([#56])
 - Fix Zipf returning values larger than `n` in rare cases ([#57])
 - Changed threshold for Ahrens-Dieter PD algorithm to >=10, matching the suggested value in the original 1982 paper
@@ -167,6 +168,7 @@ Initial release. This is equivalent to the code in `rand` 0.6.5.
 [#46]: https://github.com/rust-random/rand_distr/pull/46
 [#48]: https://github.com/rust-random/rand_distr/pull/48
 [#53]: https://github.com/rust-random/rand_distr/pull/53
+[#60]: https://github.com/rust-random/rand_distr/pull/60
 [#56]: https://github.com/rust-random/rand_distr/pull/56
 [#57]: https://github.com/rust-random/rand_distr/pull/57
 [rand#840]: https://github.com/rust-random/rand/pull/840

src/hypergeometric.rs

@@ -201,7 +201,7 @@ impl Hypergeometric {
         // J. Statist. Comput. Simul. Vol.22 (August 1985), 127-145
         // https://www.researchgate.net/publication/233212638
         const HIN_THRESHOLD: f64 = 10.0;
-        let m = ((k + 1) as f64 * (n1 + 1) as f64 / (n + 2) as f64).floor();
+        let m = ((k + 1) as f64 * (n1 + 1) as f64 / (n as f64 + 2.0)).floor();
         let sampling_method = if m - f64::max(0.0, k as f64 - n2 as f64) < HIN_THRESHOLD {
             let (initial_p, initial_x) = if k < n2 {
                 (
@@ -498,6 +498,19 @@ mod test {
         test_hypergeometric_mean_and_variance(100100, 100, 10000, &mut rng);
     }
 
+    #[test]
+    fn test_hypergeometric_near_max_population() {
+        // `n + 2` used to be computed in u64 and overflowed when the
+        // population approached u64::MAX, corrupting the mode `m` (and the
+        // sampling-method selection). Constructing and sampling must work.
+        let mut rng = crate::test::rng(738);
+        let d = Hypergeometric::new(u64::MAX - 1, 3, 2).unwrap();
+        for _ in 0..1000 {
+            let x = d.sample(&mut rng);
+            assert!(x <= 2);
+        }
+    }
+
     #[test]
     fn hypergeometric_distributions_can_be_compared() {
         assert_eq!(Hypergeometric::new(1, 2, 3), Hypergeometric::new(1, 2, 3));