Merge pull request #14531 from WarrenWeckesser/binomial-regression

BUG: random: Create a legacy implementation of random.binomial.

7 files changed, 75 insertions, 11 deletions

diff --git a/numpy/random/legacy_distributions.pxd b/numpy/random/legacy_distributions.pxd
index 7ba058054..c681388db 100644
--- a/numpy/random/legacy_distributions.pxd
+++ b/numpy/random/legacy_distributions.pxd
@@ -34,6 +34,8 @@ cdef extern from "legacy-distributions.h":
                             double nonc) nogil
     double legacy_wald(aug_bitgen_t *aug_state, double mean, double scale) nogil
     double legacy_lognormal(aug_bitgen_t *aug_state, double mean, double sigma) nogil
+    int64_t legacy_random_binomial(bitgen_t *bitgen_state, double p,
+                                   int64_t n, binomial_t *binomial) nogil
     int64_t legacy_negative_binomial(aug_bitgen_t *aug_state, double n, double p) nogil
     int64_t legacy_random_hypergeometric(bitgen_t *bitgen_state, int64_t good, int64_t bad, int64_t sample) nogil
     int64_t legacy_random_logseries(bitgen_t *bitgen_state, double p) nogil
diff --git a/numpy/random/mtrand.pyx b/numpy/random/mtrand.pyx
index 267a056e8..c469a4645 100644
--- a/numpy/random/mtrand.pyx
+++ b/numpy/random/mtrand.pyx
@@ -3086,7 +3086,9 @@ cdef class RandomState:
                 for i in range(cnt):
                     _dp = (<double*>np.PyArray_MultiIter_DATA(it, 1))[0]
                     _in = (<long*>np.PyArray_MultiIter_DATA(it, 2))[0]
-                    (<long*>np.PyArray_MultiIter_DATA(it, 0))[0] = random_binomial(&self._bitgen, _dp, _in, &self._binomial)
+                    (<long*>np.PyArray_MultiIter_DATA(it, 0))[0] = \
+                        legacy_random_binomial(&self._bitgen, _dp, _in,
+                                               &self._binomial)
 
                     np.PyArray_MultiIter_NEXT(it)
 
@@ -3099,7 +3101,8 @@ cdef class RandomState:
 
         if size is None:
             with self.lock:
-                return random_binomial(&self._bitgen, _dp, _in, &self._binomial)
+                return <long>legacy_random_binomial(&self._bitgen, _dp, _in,
+                                                    &self._binomial)
 
         randoms = <np.ndarray>np.empty(size, int)
         cnt = np.PyArray_SIZE(randoms)
@@ -3107,8 +3110,8 @@ cdef class RandomState:
 
         with self.lock, nogil:
             for i in range(cnt):
-                randoms_data[i] = random_binomial(&self._bitgen, _dp, _in,
-                                                  &self._binomial)
+                randoms_data[i] = legacy_random_binomial(&self._bitgen, _dp, _in,
+                                                         &self._binomial)
 
         return randoms
 
diff --git a/numpy/random/src/distributions/distributions.c b/numpy/random/src/distributions/distributions.c
index f8d6b9dc5..1244ffe65 100644
--- a/numpy/random/src/distributions/distributions.c
+++ b/numpy/random/src/distributions/distributions.c
@@ -901,8 +901,8 @@ RAND_INT_TYPE random_binomial_inversion(bitgen_t *bitgen_state, RAND_INT_TYPE n,
   return X;
 }
 
-RAND_INT_TYPE random_binomial(bitgen_t *bitgen_state, double p, RAND_INT_TYPE n,
-                              binomial_t *binomial) {
+int64_t random_binomial(bitgen_t *bitgen_state, double p, int64_t n,
+                        binomial_t *binomial) {
   double q;
 
   if ((n == 0LL) || (p == 0.0f))
diff --git a/numpy/random/src/distributions/distributions.h b/numpy/random/src/distributions/distributions.h
index c8cdfd20f..f2c370c07 100644
--- a/numpy/random/src/distributions/distributions.h
+++ b/numpy/random/src/distributions/distributions.h
@@ -43,11 +43,11 @@
 typedef struct s_binomial_t {
   int has_binomial; /* !=0: following parameters initialized for binomial */
   double psave;
-  int64_t nsave;
+  RAND_INT_TYPE nsave;
   double r;
   double q;
   double fm;
-  int64_t m;
+  RAND_INT_TYPE m;
   double p1;
   double xm;
   double xl;
@@ -148,8 +148,18 @@ DECLDIR double random_triangular(bitgen_t *bitgen_state, double left, double mod
 DECLDIR RAND_INT_TYPE random_poisson(bitgen_t *bitgen_state, double lam);
 DECLDIR RAND_INT_TYPE random_negative_binomial(bitgen_t *bitgen_state, double n,
                                          double p);
-DECLDIR RAND_INT_TYPE random_binomial(bitgen_t *bitgen_state, double p, RAND_INT_TYPE n,
-                                binomial_t *binomial);
+
+DECLDIR RAND_INT_TYPE random_binomial_btpe(bitgen_t *bitgen_state,
+                                           RAND_INT_TYPE n,
+                                           double p,
+                                           binomial_t *binomial);
+DECLDIR RAND_INT_TYPE random_binomial_inversion(bitgen_t *bitgen_state,
+                                                RAND_INT_TYPE n,
+                                                double p,
+                                                binomial_t *binomial);
+DECLDIR int64_t random_binomial(bitgen_t *bitgen_state, double p,
+                                int64_t n, binomial_t *binomial);
+
 DECLDIR RAND_INT_TYPE random_logseries(bitgen_t *bitgen_state, double p);
 DECLDIR RAND_INT_TYPE random_geometric_search(bitgen_t *bitgen_state, double p);
 DECLDIR RAND_INT_TYPE random_geometric_inversion(bitgen_t *bitgen_state, double p);
diff --git a/numpy/random/src/legacy/legacy-distributions.c b/numpy/random/src/legacy/legacy-distributions.c
index 4741a0352..684b3d762 100644
--- a/numpy/random/src/legacy/legacy-distributions.c
+++ b/numpy/random/src/legacy/legacy-distributions.c
@@ -215,6 +215,37 @@ double legacy_exponential(aug_bitgen_t *aug_state, double scale) {
 }
 
 
+static RAND_INT_TYPE legacy_random_binomial_original(bitgen_t *bitgen_state,
+                                                     double p,
+                                                     RAND_INT_TYPE n,
+                                                     binomial_t *binomial) {
+  double q;
+
+  if (p <= 0.5) {
+    if (p * n <= 30.0) {
+      return random_binomial_inversion(bitgen_state, n, p, binomial);
+    } else {
+      return random_binomial_btpe(bitgen_state, n, p, binomial);
+    }
+  } else {
+    q = 1.0 - p;
+    if (q * n <= 30.0) {
+      return n - random_binomial_inversion(bitgen_state, n, q, binomial);
+    } else {
+      return n - random_binomial_btpe(bitgen_state, n, q, binomial);
+    }
+  }
+}
+
+
+int64_t legacy_random_binomial(bitgen_t *bitgen_state, double p,
+                               int64_t n, binomial_t *binomial) {
+  return (int64_t) legacy_random_binomial_original(bitgen_state, p,
+                                                   (RAND_INT_TYPE) n,
+                                                   binomial);
+}
+
+
 static RAND_INT_TYPE random_hypergeometric_hyp(bitgen_t *bitgen_state,
                                                RAND_INT_TYPE good,
                                                RAND_INT_TYPE bad,
diff --git a/numpy/random/src/legacy/legacy-distributions.h b/numpy/random/src/legacy/legacy-distributions.h
index 0e58fc61c..4bc15d58e 100644
--- a/numpy/random/src/legacy/legacy-distributions.h
+++ b/numpy/random/src/legacy/legacy-distributions.h
@@ -31,6 +31,8 @@ extern double legacy_f(aug_bitgen_t *aug_state, double dfnum, double dfden);
 extern double legacy_normal(aug_bitgen_t *aug_state, double loc, double scale);
 extern double legacy_standard_gamma(aug_bitgen_t *aug_state, double shape);
 extern double legacy_exponential(aug_bitgen_t *aug_state, double scale);
+extern int64_t legacy_random_binomial(bitgen_t *bitgen_state, double p,
+                                      int64_t n, binomial_t *binomial);
 extern int64_t legacy_negative_binomial(aug_bitgen_t *aug_state, double n,
                                         double p);
 extern int64_t legacy_random_hypergeometric(bitgen_t *bitgen_state,
diff --git a/numpy/random/tests/test_randomstate_regression.py b/numpy/random/tests/test_randomstate_regression.py
index ae92319c4..edf32ea97 100644
--- a/numpy/random/tests/test_randomstate_regression.py
+++ b/numpy/random/tests/test_randomstate_regression.py
@@ -191,4 +191,20 @@ class TestRegression(object):
                              2588848963, 3684848379, 2340255427, 3638918503,
                              1819583497, 2678185683], dtype='int64')
         actual = random.randint(2**32, size=10)
-        assert_array_equal(actual, expected)
\ No newline at end of file
+        assert_array_equal(actual, expected)
+
+    def test_p_zero_stream(self):
+        # Regression test for gh-14522.  Ensure that future versions
+        # generate the same variates as version 1.16.
+        np.random.seed(12345)
+        assert_array_equal(random.binomial(1, [0, 0.25, 0.5, 0.75, 1]),
+                           [0, 0, 0, 1, 1])
+
+    def test_n_zero_stream(self):
+        # Regression test for gh-14522.  Ensure that future versions
+        # generate the same variates as version 1.16.
+        np.random.seed(8675309)
+        expected = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+                             [3, 4, 2, 3, 3, 1, 5, 3, 1, 3]])
+        assert_array_equal(random.binomial([[0], [10]], 0.25, size=(2, 10)),
+                           expected)