1 files changed, 68 insertions, 0 deletions

diff --git a/numpy/random/src/legacy/legacy-distributions.c b/numpy/random/src/legacy/legacy-distributions.c
index fd067fe8d..5b17401dd 100644
--- a/numpy/random/src/legacy/legacy-distributions.c
+++ b/numpy/random/src/legacy/legacy-distributions.c
@@ -1,3 +1,13 @@
+/*
+ * This file contains generation code for distribution that have been modified
+ * since Generator was introduced. These are preserved using identical code
+ * to what was in NumPy 1.16 so that the stream of values generated by
+ * RandomState is not changed when there are changes that affect Generator.
+ *
+ * These functions should not be changed except if they contain code that
+ * cannot be compiled. They should not be changed for bug fixes, performance
+ * improvements that can change the values produced, or enhancements to precision.
+ */
 #include "include/legacy-distributions.h"
 
 
@@ -390,3 +400,61 @@ int64_t legacy_random_logseries(bitgen_t *bitgen_state, double p) {
                                binomial_t *binomial) {
   return random_multinomial(bitgen_state, n, mnix, pix, d, binomial);
 }
+
+double legacy_vonmises(bitgen_t *bitgen_state, double mu, double kappa) {
+  double s;
+  double U, V, W, Y, Z;
+  double result, mod;
+  int neg;
+  if (npy_isnan(kappa)) {
+    return NPY_NAN;
+  }
+  if (kappa < 1e-8) {
+    return M_PI * (2 * next_double(bitgen_state) - 1);
+  } else {
+    /* with double precision rho is zero until 1.4e-8 */
+    if (kappa < 1e-5) {
+      /*
+       * second order taylor expansion around kappa = 0
+       * precise until relatively large kappas as second order is 0
+       */
+      s = (1. / kappa + kappa);
+    } else {
+        /* 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);
+    }
+
+    while (1) {
+      U = next_double(bitgen_state);
+      Z = cos(M_PI * U);
+      W = (1 + s * Z) / (s + Z);
+      Y = kappa * (s - W);
+      V = next_double(bitgen_state);
+      /*
+       * V==0.0 is ok here since Y >= 0 always leads
+       * to accept, while Y < 0 always rejects
+       */
+      if ((Y * (2 - Y) - V >= 0) || (log(Y / V) + 1 - Y >= 0)) {
+        break;
+      }
+    }
+
+    U = next_double(bitgen_state);
+
+    result = acos(W);
+    if (U < 0.5) {
+      result = -result;
+    }
+    result += mu;
+    neg = (result < 0);
+    mod = fabs(result);
+    mod = (fmod(mod + M_PI, 2 * M_PI) - M_PI);
+    if (neg) {
+      mod *= -1;
+    }
+
+    return mod;
+  }
+}
\ No newline at end of file