diff options
Diffstat (limited to 'numpy/random/src/distributions/distributions.c')
| -rw-r--r-- | numpy/random/src/distributions/distributions.c | 29 |
1 files changed, 11 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) { |