BUG: Prevent RandomState from changing

Apply vonmises fix only to Generator
Add tests for correctness

closes #17378
closes #17275

6 files changed, 45 insertions, 19 deletions

diff --git a/numpy/random/src/distributions/distributions.c b/numpy/random/src/distributions/distributions.c
index 59e11ca47..ce909933e 100644
--- a/numpy/random/src/distributions/distributions.c
+++ b/numpy/random/src/distributions/distributions.c
@@ -1,6 +1,7 @@
 #include "numpy/random/distributions.h"
 #include "ziggurat_constants.h"
 #include "logfactorial.h"
+#include <stdio.h>
 
 #if defined(_MSC_VER) && defined(_WIN64)
 #include <intrin.h>
@@ -844,8 +845,7 @@ double random_vonmises(bitgen_t *bitgen_state, double mu, double kappa) {
   }
   if (kappa < 1e-8) {
     return M_PI * (2 * next_double(bitgen_state) - 1);
-  }
-  else {
+  } else {
     /* with double precision rho is zero until 1.4e-8 */
     if (kappa < 1e-5) {
       /*
@@ -853,9 +853,12 @@ double random_vonmises(bitgen_t *bitgen_state, double mu, double kappa) {
        * precise until relatively large kappas as second order is 0
        */
       s = (1. / kappa + kappa);
-    }
-    else {
-      /* Fallback to normal distribution for big values of kappa */
+    } else {
+#ifndef NP_RANDOM_LEGACY
+      /* Fallback to normal distribution for big values of kappa
+       * Fix only applies to Generator, not RandomState
+       * RandomState defines NP_RANDOM_LEGACY=1
+       */
       if (kappa > 1e6) {
         result = mu + sqrt(1. / kappa) * random_standard_normal(bitgen_state);
         /* Check if result is within bounds */
@@ -866,12 +869,15 @@ double random_vonmises(bitgen_t *bitgen_state, double mu, double kappa) {
           return result - 2*M_PI;
         }
         return result;
-      }
-      else {
+      } else {
+#endif
+        /* Path for 1e-5 <= kappa <= 1e6 */
         double r = 1 + sqrt(1 + 4 * kappa * kappa);
         double rho = (r - sqrt(2 * r)) / (2 * kappa);
         s = (1 + rho * rho) / (2 * rho);
+#if !defined NP_RANDOM_LEGACY
       }
+#endif
     }
 
     while (1) {
diff --git a/numpy/random/tests/test_generator_mt19937.py b/numpy/random/tests/test_generator_mt19937.py
index d68bcd38b..f6bd985b4 100644
--- a/numpy/random/tests/test_generator_mt19937.py
+++ b/numpy/random/tests/test_generator_mt19937.py
@@ -10,7 +10,7 @@ from numpy.testing import (
     assert_warns, assert_no_warnings, assert_array_equal,
     assert_array_almost_equal, suppress_warnings)
 
-from numpy.random import Generator, MT19937, SeedSequence
+from numpy.random import Generator, MT19937, SeedSequence, RandomState
 
 random = Generator(MT19937())
 
@@ -1735,6 +1735,26 @@ class TestRandomDist:
         r = random.vonmises(mu=0., kappa=np.nan)
         assert_(np.isnan(r))
 
+    @pytest.mark.parametrize("kappa", [1e4, 1e15])
+    def test_vonmises_large_kappa(self, kappa):
+        random = Generator(MT19937(self.seed))
+        rs = RandomState(random.bit_generator)
+        state = random.bit_generator.state
+
+        random_state_vals = rs.vonmises(0, kappa, size=10)
+        random.bit_generator.state = state
+        gen_vals = random.vonmises(0, kappa, size=10)
+        if kappa < 1e6:
+            assert_allclose(random_state_vals, gen_vals)
+        else:
+            assert np.all(random_state_vals != gen_vals)
+
+    @pytest.mark.parametrize("mu", [-7., -np.pi, -3.1, np.pi, 3.2])
+    @pytest.mark.parametrize("kappa", [1e-9, 1e-6, 1, 1e3, 1e15])
+    def test_vonmises_large_kappa_range(self, mu, kappa):
+        r = random.vonmises(mu, kappa, 50)
+        assert_(np.all(r > -np.pi) and np.all(r <= np.pi))
+
     def test_wald(self):
         random = Generator(MT19937(self.seed))
         actual = random.wald(mean=1.23, scale=1.54, size=(3, 2))
diff --git a/numpy/random/tests/test_generator_mt19937_regressions.py b/numpy/random/tests/test_generator_mt19937_regressions.py
index 314713cbe..9f6dcdc6b 100644
--- a/numpy/random/tests/test_generator_mt19937_regressions.py
+++ b/numpy/random/tests/test_generator_mt19937_regressions.py
@@ -1,22 +1,19 @@
 from numpy.testing import (assert_, assert_array_equal)
 import numpy as np
 import pytest
-from numpy.random import Generator, MT19937
+from numpy.random import Generator, MT19937, RandomState
 
 mt19937 = Generator(MT19937())
 
 
 class TestRegression:
 
-    def test_VonMises_range(self):
+    def test_vonmises_range(self):
         # Make sure generated random variables are in [-pi, pi].
         # Regression test for ticket #986.
         for mu in np.linspace(-7., 7., 5):
             r = mt19937.vonmises(mu, 1, 50)
             assert_(np.all(r > -np.pi) and np.all(r <= np.pi))
-        for mu in [-np.pi, np.pi]:
-            r = mt19937.vonmises(mu, 10**6 + 1, 50)
-            assert_(np.all(r > -np.pi) and np.all(r <= np.pi))
 
     def test_hypergeometric_range(self):
         # Test for ticket #921
diff --git a/numpy/random/tests/test_randomstate.py b/numpy/random/tests/test_randomstate.py
index 7f5f08050..b16275b70 100644
--- a/numpy/random/tests/test_randomstate.py
+++ b/numpy/random/tests/test_randomstate.py
@@ -1238,6 +1238,15 @@ class TestRandomDist:
         r = random.vonmises(mu=0., kappa=1.1e-8, size=10**6)
         assert_(np.isfinite(r).all())
 
+    def test_vonmises_large(self):
+        # guard against changes in RandomState when Generator is fixed
+        random.seed(self.seed)
+        actual = random.vonmises(mu=0., kappa=1e7, size=3)
+        desired = np.array([4.634253748521111e-04,
+                            3.558873596114509e-04,
+                            -2.337119622577433e-04])
+        assert_array_almost_equal(actual, desired, decimal=8)
+
     def test_vonmises_nan(self):
         random.seed(self.seed)
         r = random.vonmises(mu=0., kappa=np.nan)
diff --git a/numpy/random/tests/test_randomstate_regression.py b/numpy/random/tests/test_randomstate_regression.py
index d7ca540d8..0bf361e5e 100644
--- a/numpy/random/tests/test_randomstate_regression.py
+++ b/numpy/random/tests/test_randomstate_regression.py
@@ -18,9 +18,6 @@ class TestRegression:
         for mu in np.linspace(-7., 7., 5):
             r = random.vonmises(mu, 1, 50)
             assert_(np.all(r > -np.pi) and np.all(r <= np.pi))
-        for mu in [-np.pi, np.pi]:
-            r = random.vonmises(mu, 10**6 + 1, 50)
-            assert_(np.all(r > -np.pi) and np.all(r <= np.pi))
 
     def test_hypergeometric_range(self):
         # Test for ticket #921
diff --git a/numpy/random/tests/test_regression.py b/numpy/random/tests/test_regression.py
index 2ebd4c12f..54d5a3efb 100644
--- a/numpy/random/tests/test_regression.py
+++ b/numpy/random/tests/test_regression.py
@@ -14,9 +14,6 @@ class TestRegression:
         for mu in np.linspace(-7., 7., 5):
             r = random.mtrand.vonmises(mu, 1, 50)
             assert_(np.all(r > -np.pi) and np.all(r <= np.pi))
-        for mu in [-np.pi, np.pi]:
-            r = random.mtrand.vonmises(mu, 10**6 + 1, 50)
-            assert_(np.all(r > -np.pi) and np.all(r <= np.pi))
 
     def test_hypergeometric_range(self):
         # Test for ticket #921