Merge pull request #18666 from bashtage/rayleigh-perf

ENH: Use exponentials in place of inversion in Rayleigh and geometric

6 files changed, 29 insertions, 10 deletions

diff --git a/doc/release/upcoming_changes/18666.improvement.rst b/doc/release/upcoming_changes/18666.improvement.rst
new file mode 100644
index 000000000..70b87ecf4
--- /dev/null
+++ b/doc/release/upcoming_changes/18666.improvement.rst
@@ -0,0 +1,9 @@
+``Generator.rayleigh`` and ``Generator.geometric`` performance improved
+-----------------------------------------------------------------------
+The performance of Rayleigh and geometric random variate generation
+in ``Generator`` has improved. These are both transformation of exponential
+random variables and the slow log-based inverse cdf transformation has
+been replaced with the Ziggurat-based exponential variate generator.
+
+This change breaks the stream of variates generated  when variates from
+either of these distributions are produced.
diff --git a/numpy/random/include/legacy-distributions.h b/numpy/random/include/legacy-distributions.h
index 3d882b73b..f6c5cf053 100644
--- a/numpy/random/include/legacy-distributions.h
+++ b/numpy/random/include/legacy-distributions.h
@@ -17,6 +17,7 @@ extern double legacy_weibull(aug_bitgen_t *aug_state, double a);
 extern double legacy_power(aug_bitgen_t *aug_state, double a);
 extern double legacy_gamma(aug_bitgen_t *aug_state, double shape, double scale);
 extern double legacy_chisquare(aug_bitgen_t *aug_state, double df);
+extern double legacy_rayleigh(bitgen_t *bitgen_state, double mode);
 extern double legacy_noncentral_chisquare(aug_bitgen_t *aug_state, double df,
                                           double nonc);
 extern double legacy_noncentral_f(aug_bitgen_t *aug_state, double dfnum,
diff --git a/numpy/random/mtrand.pyx b/numpy/random/mtrand.pyx
index 23cb5ea31..e166634be 100644
--- a/numpy/random/mtrand.pyx
+++ b/numpy/random/mtrand.pyx
@@ -79,6 +79,7 @@ cdef extern from "include/legacy-distributions.h":
     double legacy_gamma(aug_bitgen_t *aug_state, double shape, double scale) nogil
     double legacy_power(aug_bitgen_t *aug_state, double a) nogil
     double legacy_chisquare(aug_bitgen_t *aug_state, double df) nogil
+    double legacy_rayleigh(aug_bitgen_t *aug_state, double mode) nogil
     double legacy_noncentral_chisquare(aug_bitgen_t *aug_state, double df,
                                     double nonc) nogil
     double legacy_noncentral_f(aug_bitgen_t *aug_state, double dfnum, double dfden,
@@ -3086,7 +3087,7 @@ cdef class RandomState:
         0.087300000000000003 # random
 
         """
-        return cont(&random_rayleigh, &self._bitgen, size, self.lock, 1,
+        return cont(&legacy_rayleigh, &self._bitgen, size, self.lock, 1,
                     scale, 'scale', CONS_NON_NEGATIVE,
                     0.0, '', CONS_NONE,
                     0.0, '', CONS_NONE, None)
diff --git a/numpy/random/src/distributions/distributions.c b/numpy/random/src/distributions/distributions.c
index 6b4deb925..9bdfa9bea 100644
--- a/numpy/random/src/distributions/distributions.c
+++ b/numpy/random/src/distributions/distributions.c
@@ -508,7 +508,7 @@ double random_lognormal(bitgen_t *bitgen_state, double mean, double sigma) {
 }
 
 double random_rayleigh(bitgen_t *bitgen_state, double mode) {
-  return mode * sqrt(-2.0 * npy_log1p(-next_double(bitgen_state)));
+  return mode * sqrt(2.0 * random_standard_exponential(bitgen_state));
 }
 
 double random_standard_t(bitgen_t *bitgen_state, double df) {
@@ -960,7 +960,7 @@ RAND_INT_TYPE random_geometric_search(bitgen_t *bitgen_state, double p) {
 }
 
 int64_t random_geometric_inversion(bitgen_t *bitgen_state, double p) {
-  return (int64_t)ceil(npy_log1p(-next_double(bitgen_state)) / npy_log1p(-p));
+  return (int64_t)ceil(-random_standard_exponential(bitgen_state) / npy_log1p(-p));
 }
 
 int64_t random_geometric(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 bfea15e40..443c1a4bf 100644
--- a/numpy/random/src/legacy/legacy-distributions.c
+++ b/numpy/random/src/legacy/legacy-distributions.c
@@ -112,6 +112,10 @@ double legacy_chisquare(aug_bitgen_t *aug_state, double df) {
   return 2.0 * legacy_standard_gamma(aug_state, df / 2.0);
 }
 
+double legacy_rayleigh(bitgen_t *bitgen_state, double mode) {
+  return mode * sqrt(-2.0 * npy_log1p(-next_double(bitgen_state)));
+}
+
 double legacy_noncentral_chisquare(aug_bitgen_t *aug_state, double df,
                                    double nonc) {
   double out;
diff --git a/numpy/random/tests/test_generator_mt19937.py b/numpy/random/tests/test_generator_mt19937.py
index 310545e0d..0108d84b3 100644
--- a/numpy/random/tests/test_generator_mt19937.py
+++ b/numpy/random/tests/test_generator_mt19937.py
@@ -1252,9 +1252,9 @@ class TestRandomDist:
     def test_geometric(self):
         random = Generator(MT19937(self.seed))
         actual = random.geometric(.123456789, size=(3, 2))
-        desired = np.array([[ 1, 10],
-                            [ 1, 12],
-                            [ 9, 10]])
+        desired = np.array([[1, 11],
+                            [1, 12],
+                            [11, 17]])
         assert_array_equal(actual, desired)
 
     def test_geometric_exceptions(self):
@@ -1557,9 +1557,9 @@ class TestRandomDist:
     def test_rayleigh(self):
         random = Generator(MT19937(self.seed))
         actual = random.rayleigh(scale=10, size=(3, 2))
-        desired = np.array([[ 4.51734079831581, 15.6802442485758 ],
-                            [ 4.19850651287094, 17.08718809823704],
-                            [14.7907457708776 , 15.85545333419775]])
+        desired = np.array([[4.19494429102666, 16.66920198906598],
+                            [3.67184544902662, 17.74695521962917],
+                            [16.27935397855501, 21.08355560691792]])
         assert_array_almost_equal(actual, desired, decimal=14)
 
     def test_rayleigh_0(self):
@@ -2114,7 +2114,11 @@ class TestBroadcast:
     def test_rayleigh(self):
         scale = [1]
         bad_scale = [-1]
-        desired = np.array([0.60439534475066, 0.66120048396359, 1.67873398389499])
+        desired = np.array(
+            [1.1597068009872629,
+             0.6539188836253857,
+             1.1981526554349398]
+        )
 
         random = Generator(MT19937(self.seed))
         actual = random.rayleigh(scale * 3)