BUG: Protect gamma generation from 0 input

Add protection for 0 input for gamma random variables
Add protection in legacy for gamma random variables
Add protection for 0 input in Weibull
Add test of 0 protection
Ensure tests run on legacy generators for bad values are also run on the
current set of generators
Use same variable definition in declaration and function
Sync doc changes to Wald and Noncentral Chi2

7 files changed, 72 insertions, 14 deletions

diff --git a/_randomgen/randomgen/generator.pyx b/_randomgen/randomgen/generator.pyx
index 2d45f8aca..a4e3dd3a1 100644
--- a/_randomgen/randomgen/generator.pyx
+++ b/_randomgen/randomgen/generator.pyx
@@ -1954,17 +1954,9 @@ cdef class RandomGenerator:
 
         where :math:`Y_{q}` is the Chi-square with q degrees of freedom.
 
-        In Delhi (2007), it is noted that the noncentral chi-square is
-        useful in bombing and coverage problems, the probability of
-        killing the point target given by the noncentral chi-squared
-        distribution.
-
         References
         ----------
-        .. [1] Delhi, M.S. Holla, "On a noncentral chi-square distribution in
-               the analysis of weapon systems effectiveness", Metrika,
-               Volume 15, Number 1 / December, 1970.
-        .. [2] Wikipedia, "Noncentral chi-square distribution"
+        .. [1] Wikipedia, "Noncentral chi-square distribution"
                https://en.wikipedia.org/wiki/Noncentral_chi-square_distribution
 
         Examples
@@ -3016,9 +3008,9 @@ cdef class RandomGenerator:
         Parameters
         ----------
         mean : float or array_like of floats
-            Distribution mean, should be > 0.
+            Distribution mean, must be > 0.
         scale : float or array_like of floats
-            Scale parameter, should be >= 0.
+            Scale parameter, must be > 0.
         size : int or tuple of ints, optional
             Output shape.  If the given shape is, e.g., ``(m, n, k)``, then
             ``m * n * k`` samples are drawn.  If size is ``None`` (default),
diff --git a/_randomgen/randomgen/src/distributions/distributions.c b/_randomgen/randomgen/src/distributions/distributions.c
index 4e7493afa..a9d8a308d 100644
--- a/_randomgen/randomgen/src/distributions/distributions.c
+++ b/_randomgen/randomgen/src/distributions/distributions.c
@@ -346,6 +346,8 @@ static NPY_INLINE double standard_gamma_zig(brng_t *brng_state, double shape) {
 
   if (shape == 1.0) {
     return random_standard_exponential_zig(brng_state);
+  } else if (shape == 0.0) {
+    return 0.0;
   } else if (shape < 1.0) {
     for (;;) {
       U = next_double(brng_state);
@@ -388,6 +390,8 @@ static NPY_INLINE float standard_gamma_zig_f(brng_t *brng_state, float shape) {
 
   if (shape == 1.0f) {
     return random_standard_exponential_zig_f(brng_state);
+  } else if (shape == 0.0) {
+    return 0.0;
   } else if (shape < 1.0f) {
     for (;;) {
       U = next_float(brng_state);
diff --git a/_randomgen/randomgen/src/distributions/distributions.h b/_randomgen/randomgen/src/distributions/distributions.h
index e618052c1..5cf9c72b2 100644
--- a/_randomgen/randomgen/src/distributions/distributions.h
+++ b/_randomgen/randomgen/src/distributions/distributions.h
@@ -210,7 +210,7 @@ DECLDIR void random_bounded_uint32_fill(brng_t *brng_state, uint32_t off,
 DECLDIR void random_bounded_uint16_fill(brng_t *brng_state, uint16_t off,
                                         uint16_t rng, npy_intp cnt,
                                         bool use_masked, uint16_t *out);
-DECLDIR void random_bounded_uint8_fill(brng_t *brng_state, npy_bool off,
+DECLDIR void random_bounded_uint8_fill(brng_t *brng_state, uint8_t off,
                                        uint8_t rng, npy_intp cnt,
                                        bool use_masked, uint8_t *out);
 DECLDIR void random_bounded_bool_fill(brng_t *brng_state, npy_bool off,
diff --git a/_randomgen/randomgen/src/legacy/distributions-boxmuller.c b/_randomgen/randomgen/src/legacy/distributions-boxmuller.c
index 7ad1742cc..768de066c 100644
--- a/_randomgen/randomgen/src/legacy/distributions-boxmuller.c
+++ b/_randomgen/randomgen/src/legacy/distributions-boxmuller.c
@@ -39,6 +39,9 @@ double legacy_standard_gamma(aug_brng_t *aug_state, double shape) {
 
   if (shape == 1.0) {
     return legacy_standard_exponential(aug_state);
+  }
+  else if (shape == 0.0) {
+    return 0.0;
   } else if (shape < 1.0) {
     for (;;) {
       U = legacy_double(aug_state);
@@ -84,6 +87,9 @@ double legacy_pareto(aug_brng_t *aug_state, double a) {
 }
 
 double legacy_weibull(aug_brng_t *aug_state, double a) {
+  if (a == 0.0) {
+    return 0.0;
+  }
   return pow(legacy_standard_exponential(aug_state), 1. / a);
 }
 
diff --git a/_randomgen/randomgen/tests/test_legacy.py b/_randomgen/randomgen/tests/test_legacy.py
index e45ba3619..21c56946f 100644
--- a/_randomgen/randomgen/tests/test_legacy.py
+++ b/_randomgen/randomgen/tests/test_legacy.py
@@ -10,3 +10,8 @@ def test_pickle():
     lg2 = pickle.loads(pickle.dumps(lg))
     assert lg.standard_normal() == lg2.standard_normal()
     assert lg.random_sample() == lg2.random_sample()
+
+
+def test_weibull():
+    lg = LegacyGenerator()
+    assert lg.weibull(0.0) == 0.0
diff --git a/_randomgen/randomgen/tests/test_numpy_mt19937.py b/_randomgen/randomgen/tests/test_numpy_mt19937.py
index a43d7d601..df13820ca 100644
--- a/_randomgen/randomgen/tests/test_numpy_mt19937.py
+++ b/_randomgen/randomgen/tests/test_numpy_mt19937.py
@@ -1161,11 +1161,13 @@ class TestBroadcast(object):
         actual = normal(loc * 3, scale)
         assert_array_almost_equal(actual, desired, decimal=14)
         assert_raises(ValueError, normal, loc * 3, bad_scale)
+        assert_raises(ValueError, mt19937.normal, loc * 3, bad_scale)
 
         self.set_seed()
         actual = normal(loc, scale * 3)
         assert_array_almost_equal(actual, desired, decimal=14)
         assert_raises(ValueError, normal, loc, bad_scale * 3)
+        assert_raises(ValueError, mt19937.normal, loc, bad_scale * 3)
 
     def test_beta(self):
         a = [1]
@@ -1182,12 +1184,14 @@ class TestBroadcast(object):
         assert_array_almost_equal(actual, desired, decimal=14)
         assert_raises(ValueError, beta, bad_a * 3, b)
         assert_raises(ValueError, beta, a * 3, bad_b)
+        assert_raises(ValueError, mt19937.beta, bad_a * 3, b)
+        assert_raises(ValueError, mt19937.beta, a * 3, bad_b)
 
         self.set_seed()
         actual = beta(a, b * 3)
         assert_array_almost_equal(actual, desired, decimal=14)
-        assert_raises(ValueError, beta, bad_a, b * 3)
-        assert_raises(ValueError, beta, a, bad_b * 3)
+        assert_raises(ValueError, mt19937.beta, bad_a, b * 3)
+        assert_raises(ValueError, mt19937.beta, a, bad_b * 3)
 
     def test_exponential(self):
         scale = [1]
@@ -1201,6 +1205,7 @@ class TestBroadcast(object):
         actual = exponential(scale * 3)
         assert_array_almost_equal(actual, desired, decimal=14)
         assert_raises(ValueError, exponential, bad_scale * 3)
+        assert_raises(ValueError, mt19937.exponential, bad_scale * 3)
 
     def test_standard_gamma(self):
         shape = [1]
@@ -1214,6 +1219,7 @@ class TestBroadcast(object):
         actual = std_gamma(shape * 3)
         assert_array_almost_equal(actual, desired, decimal=14)
         assert_raises(ValueError, std_gamma, bad_shape * 3)
+        assert_raises(ValueError, mt19937.standard_gamma, bad_shape * 3)
 
     def test_gamma(self):
         shape = [1]
@@ -1230,12 +1236,16 @@ class TestBroadcast(object):
         assert_array_almost_equal(actual, desired, decimal=14)
         assert_raises(ValueError, gamma, bad_shape * 3, scale)
         assert_raises(ValueError, gamma, shape * 3, bad_scale)
+        assert_raises(ValueError, mt19937.gamma, bad_shape * 3, scale)
+        assert_raises(ValueError, mt19937.gamma, shape * 3, bad_scale)
 
         self.set_seed()
         actual = gamma(shape, scale * 3)
         assert_array_almost_equal(actual, desired, decimal=14)
         assert_raises(ValueError, gamma, bad_shape, scale * 3)
         assert_raises(ValueError, gamma, shape, bad_scale * 3)
+        assert_raises(ValueError, mt19937.gamma, bad_shape, scale * 3)
+        assert_raises(ValueError, mt19937.gamma, shape, bad_scale * 3)
 
     def test_f(self):
         dfnum = [1]
@@ -1252,12 +1262,16 @@ class TestBroadcast(object):
         assert_array_almost_equal(actual, desired, decimal=14)
         assert_raises(ValueError, f, bad_dfnum * 3, dfden)
         assert_raises(ValueError, f, dfnum * 3, bad_dfden)
+        assert_raises(ValueError, mt19937.f, bad_dfnum * 3, dfden)
+        assert_raises(ValueError, mt19937.f, dfnum * 3, bad_dfden)
 
         self.set_seed()
         actual = f(dfnum, dfden * 3)
         assert_array_almost_equal(actual, desired, decimal=14)
         assert_raises(ValueError, f, bad_dfnum, dfden * 3)
         assert_raises(ValueError, f, dfnum, bad_dfden * 3)
+        assert_raises(ValueError, mt19937.f, bad_dfnum, dfden * 3)
+        assert_raises(ValueError, mt19937.f, dfnum, bad_dfden * 3)
 
     def test_noncentral_f(self):
         dfnum = [2]
@@ -1277,6 +1291,9 @@ class TestBroadcast(object):
         assert_raises(ValueError, nonc_f, bad_dfnum * 3, dfden, nonc)
         assert_raises(ValueError, nonc_f, dfnum * 3, bad_dfden, nonc)
         assert_raises(ValueError, nonc_f, dfnum * 3, dfden, bad_nonc)
+        assert_raises(ValueError, mt19937.noncentral_f, bad_dfnum * 3, dfden, nonc)
+        assert_raises(ValueError, mt19937.noncentral_f, dfnum * 3, bad_dfden, nonc)
+        assert_raises(ValueError, mt19937.noncentral_f, dfnum * 3, dfden, bad_nonc)
 
         self.set_seed()
         actual = nonc_f(dfnum, dfden * 3, nonc)
@@ -1284,6 +1301,9 @@ class TestBroadcast(object):
         assert_raises(ValueError, nonc_f, bad_dfnum, dfden * 3, nonc)
         assert_raises(ValueError, nonc_f, dfnum, bad_dfden * 3, nonc)
         assert_raises(ValueError, nonc_f, dfnum, dfden * 3, bad_nonc)
+        assert_raises(ValueError, mt19937.noncentral_f, bad_dfnum, dfden * 3, nonc)
+        assert_raises(ValueError, mt19937.noncentral_f, dfnum, bad_dfden * 3, nonc)
+        assert_raises(ValueError, mt19937.noncentral_f, dfnum, dfden * 3, bad_nonc)
 
         self.set_seed()
         actual = nonc_f(dfnum, dfden, nonc * 3)
@@ -1291,6 +1311,9 @@ class TestBroadcast(object):
         assert_raises(ValueError, nonc_f, bad_dfnum, dfden, nonc * 3)
         assert_raises(ValueError, nonc_f, dfnum, bad_dfden, nonc * 3)
         assert_raises(ValueError, nonc_f, dfnum, dfden, bad_nonc * 3)
+        assert_raises(ValueError, mt19937.noncentral_f, bad_dfnum, dfden, nonc * 3)
+        assert_raises(ValueError, mt19937.noncentral_f, dfnum, bad_dfden, nonc * 3)
+        assert_raises(ValueError, mt19937.noncentral_f, dfnum, dfden, bad_nonc * 3)
 
     def test_chisquare(self):
         df = [1]
@@ -1320,12 +1343,16 @@ class TestBroadcast(object):
         assert_array_almost_equal(actual, desired, decimal=14)
         assert_raises(ValueError, nonc_chi, bad_df * 3, nonc)
         assert_raises(ValueError, nonc_chi, df * 3, bad_nonc)
+        assert_raises(ValueError, mt19937.noncentral_chisquare, bad_df * 3, nonc)
+        assert_raises(ValueError, mt19937.noncentral_chisquare, df * 3, bad_nonc)
 
         self.set_seed()
         actual = nonc_chi(df, nonc * 3)
         assert_array_almost_equal(actual, desired, decimal=14)
         assert_raises(ValueError, nonc_chi, bad_df, nonc * 3)
         assert_raises(ValueError, nonc_chi, df, bad_nonc * 3)
+        assert_raises(ValueError, mt19937.noncentral_chisquare, bad_df, nonc * 3)
+        assert_raises(ValueError, mt19937.noncentral_chisquare, df, bad_nonc * 3)
 
     def test_standard_t(self):
         df = [1]
@@ -1339,6 +1366,7 @@ class TestBroadcast(object):
         actual = t(df * 3)
         assert_array_almost_equal(actual, desired, decimal=14)
         assert_raises(ValueError, t, bad_df * 3)
+        assert_raises(ValueError, mt19937.standard_t, bad_df * 3)
 
     def test_vonmises(self):
         mu = [2]
@@ -1371,6 +1399,7 @@ class TestBroadcast(object):
         actual = pareto(a * 3)
         assert_array_almost_equal(actual, desired, decimal=14)
         assert_raises(ValueError, pareto, bad_a * 3)
+        assert_raises(ValueError, mt19937.pareto, bad_a * 3)
 
     def test_weibull(self):
         a = [1]
@@ -1384,6 +1413,7 @@ class TestBroadcast(object):
         actual = weibull(a * 3)
         assert_array_almost_equal(actual, desired, decimal=14)
         assert_raises(ValueError, weibull, bad_a * 3)
+        assert_raises(ValueError, mt19937.weibull, bad_a * 3)
 
     def test_power(self):
         a = [1]
@@ -1397,6 +1427,7 @@ class TestBroadcast(object):
         actual = power(a * 3)
         assert_array_almost_equal(actual, desired, decimal=14)
         assert_raises(ValueError, power, bad_a * 3)
+        assert_raises(ValueError, mt19937.power, bad_a * 3)
 
     def test_laplace(self):
         loc = [0]
@@ -1468,11 +1499,13 @@ class TestBroadcast(object):
         actual = lognormal(mean * 3, sigma)
         assert_array_almost_equal(actual, desired, decimal=14)
         assert_raises(ValueError, lognormal, mean * 3, bad_sigma)
+        assert_raises(ValueError, mt19937.lognormal, mean * 3, bad_sigma)
 
         self.set_seed()
         actual = lognormal(mean, sigma * 3)
         assert_array_almost_equal(actual, desired, decimal=14)
         assert_raises(ValueError, lognormal, mean, bad_sigma * 3)
+        assert_raises(ValueError, mt19937.lognormal, mean, bad_sigma * 3)
 
     def test_rayleigh(self):
         scale = [1]
@@ -1502,12 +1535,16 @@ class TestBroadcast(object):
         assert_array_almost_equal(actual, desired, decimal=14)
         assert_raises(ValueError, wald, bad_mean * 3, scale)
         assert_raises(ValueError, wald, mean * 3, bad_scale)
+        assert_raises(ValueError, mt19937.wald, bad_mean * 3, scale)
+        assert_raises(ValueError, mt19937.wald, mean * 3, bad_scale)
 
         self.set_seed()
         actual = wald(mean, scale * 3)
         assert_array_almost_equal(actual, desired, decimal=14)
         assert_raises(ValueError, wald, bad_mean, scale * 3)
         assert_raises(ValueError, wald, mean, bad_scale * 3)
+        assert_raises(ValueError, mt19937.wald, bad_mean, scale * 3)
+        assert_raises(ValueError, mt19937.wald, mean, bad_scale * 3)
 
     def test_triangular(self):
         left = [1]
@@ -1583,6 +1620,9 @@ class TestBroadcast(object):
         assert_raises(ValueError, neg_binom, bad_n * 3, p)
         assert_raises(ValueError, neg_binom, n * 3, bad_p_one)
         assert_raises(ValueError, neg_binom, n * 3, bad_p_two)
+        assert_raises(ValueError, mt19937.negative_binomial, bad_n * 3, p)
+        assert_raises(ValueError, mt19937.negative_binomial, n * 3, bad_p_one)
+        assert_raises(ValueError, mt19937.negative_binomial, n * 3, bad_p_two)
 
         self.set_seed()
         actual = neg_binom(n, p * 3)
@@ -1590,6 +1630,9 @@ class TestBroadcast(object):
         assert_raises(ValueError, neg_binom, bad_n, p * 3)
         assert_raises(ValueError, neg_binom, n, bad_p_one * 3)
         assert_raises(ValueError, neg_binom, n, bad_p_two * 3)
+        assert_raises(ValueError, mt19937.negative_binomial, bad_n, p * 3)
+        assert_raises(ValueError, mt19937.negative_binomial, n, bad_p_one * 3)
+        assert_raises(ValueError, mt19937.negative_binomial, n, bad_p_two * 3)
 
     def test_poisson(self):
         max_lam = random.poisson_lam_max
diff --git a/_randomgen/randomgen/tests/test_numpy_mt19937_regressions.py b/_randomgen/randomgen/tests/test_numpy_mt19937_regressions.py
index 4e51327aa..bc644e122 100644
--- a/_randomgen/randomgen/tests/test_numpy_mt19937_regressions.py
+++ b/_randomgen/randomgen/tests/test_numpy_mt19937_regressions.py
@@ -158,3 +158,11 @@ class TestRegression(object):
         perm = mt19937.permutation(m)
         assert_array_equal(perm, np.array([2, 1, 4, 0, 3]))
         assert_array_equal(m.__array__(), np.arange(5))
+
+    def test_gamma_0(self):
+        assert mt19937.standard_gamma(0.0) == 0.0
+        assert_array_equal(mt19937.standard_gamma([0.0]), 0.0)
+
+        actual = mt19937.standard_gamma([0.0], dtype='float')
+        expected = np.array([0.], dtype=np.float32)
+        assert_array_equal(actual, expected)