CLN: Move to legacy function

Avoid conditional compilation and move old version to legacy_vonmises
Small clean up
Additional comments

2 files changed, 79 insertions, 18 deletions

diff --git a/numpy/random/src/distributions/distributions.c b/numpy/random/src/distributions/distributions.c
index ce909933e..f47c54a53 100644
--- a/numpy/random/src/distributions/distributions.c
+++ b/numpy/random/src/distributions/distributions.c
@@ -1,7 +1,6 @@
 #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,6 +843,7 @@ double random_vonmises(bitgen_t *bitgen_state, double mu, double kappa) {
     return NPY_NAN;
   }
   if (kappa < 1e-8) {
+    /* Use a uniform for very small values of kappa */
     return M_PI * (2 * next_double(bitgen_state) - 1);
   } else {
     /* with double precision rho is zero until 1.4e-8 */
@@ -854,30 +854,23 @@ double random_vonmises(bitgen_t *bitgen_state, double mu, double kappa) {
        */
       s = (1. / kappa + 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) {
+      if (kappa <= 1e6) {
+        /* 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);
+      } else {
+        /* Fallback to wrapped normal distribution for kappa > 1e6 */
         result = mu + sqrt(1. / kappa) * random_standard_normal(bitgen_state);
-        /* Check if result is within bounds */
+        /* Ensure result is within bounds */
         if (result < -M_PI) {
-          return result + 2*M_PI;
+          result += 2*M_PI;
         }
         if (result > M_PI) {
-          return result - 2*M_PI;
+          result -= 2*M_PI;
         }
         return result;
-      } 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/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