diff options
Diffstat (limited to 'numpy/random/src')
| -rw-r--r-- | numpy/random/src/aligned_malloc/aligned_malloc.c | 9 | ||||
| -rw-r--r-- | numpy/random/src/aligned_malloc/aligned_malloc.h | 54 | ||||
| -rw-r--r-- | numpy/random/src/distributions/distributions.c | 303 | ||||
| -rw-r--r-- | numpy/random/src/distributions/distributions.h | 205 | ||||
| -rw-r--r-- | numpy/random/src/distributions/random_hypergeometric.c | 8 | ||||
| -rw-r--r-- | numpy/random/src/distributions/random_mvhg_count.c | 131 | ||||
| -rw-r--r-- | numpy/random/src/distributions/random_mvhg_marginals.c | 138 | ||||
| -rw-r--r-- | numpy/random/src/entropy/entropy.c | 114 | ||||
| -rw-r--r-- | numpy/random/src/entropy/entropy.h | 14 | ||||
| -rw-r--r-- | numpy/random/src/legacy/legacy-distributions.c | 41 | ||||
| -rw-r--r-- | numpy/random/src/legacy/legacy-distributions.h | 52 | ||||
| -rw-r--r-- | numpy/random/src/philox/philox.h | 2 | ||||
| -rw-r--r-- | numpy/random/src/sfc64/sfc64.h | 2 |
13 files changed, 399 insertions, 674 deletions
diff --git a/numpy/random/src/aligned_malloc/aligned_malloc.c b/numpy/random/src/aligned_malloc/aligned_malloc.c deleted file mode 100644 index 6e8192cfb..000000000 --- a/numpy/random/src/aligned_malloc/aligned_malloc.c +++ /dev/null @@ -1,9 +0,0 @@ -#include "aligned_malloc.h" - -static NPY_INLINE void *PyArray_realloc_aligned(void *p, size_t n); - -static NPY_INLINE void *PyArray_malloc_aligned(size_t n); - -static NPY_INLINE void *PyArray_calloc_aligned(size_t n, size_t s); - -static NPY_INLINE void PyArray_free_aligned(void *p);
\ No newline at end of file diff --git a/numpy/random/src/aligned_malloc/aligned_malloc.h b/numpy/random/src/aligned_malloc/aligned_malloc.h deleted file mode 100644 index ea24f6d23..000000000 --- a/numpy/random/src/aligned_malloc/aligned_malloc.h +++ /dev/null @@ -1,54 +0,0 @@ -#ifndef _RANDOMDGEN__ALIGNED_MALLOC_H_ -#define _RANDOMDGEN__ALIGNED_MALLOC_H_ - -#include "Python.h" -#include "numpy/npy_common.h" - -#define NPY_MEMALIGN 16 /* 16 for SSE2, 32 for AVX, 64 for Xeon Phi */ - -static NPY_INLINE void *PyArray_realloc_aligned(void *p, size_t n) -{ - void *p1, **p2, *base; - size_t old_offs, offs = NPY_MEMALIGN - 1 + sizeof(void *); - if (NPY_UNLIKELY(p != NULL)) - { - base = *(((void **)p) - 1); - if (NPY_UNLIKELY((p1 = PyMem_Realloc(base, n + offs)) == NULL)) - return NULL; - if (NPY_LIKELY(p1 == base)) - return p; - p2 = (void **)(((Py_uintptr_t)(p1) + offs) & ~(NPY_MEMALIGN - 1)); - old_offs = (size_t)((Py_uintptr_t)p - (Py_uintptr_t)base); - memmove((void *)p2, ((char *)p1) + old_offs, n); - } - else - { - if (NPY_UNLIKELY((p1 = PyMem_Malloc(n + offs)) == NULL)) - return NULL; - p2 = (void **)(((Py_uintptr_t)(p1) + offs) & ~(NPY_MEMALIGN - 1)); - } - *(p2 - 1) = p1; - return (void *)p2; -} - -static NPY_INLINE void *PyArray_malloc_aligned(size_t n) -{ - return PyArray_realloc_aligned(NULL, n); -} - -static NPY_INLINE void *PyArray_calloc_aligned(size_t n, size_t s) -{ - void *p; - if (NPY_UNLIKELY((p = PyArray_realloc_aligned(NULL, n * s)) == NULL)) - return NULL; - memset(p, 0, n * s); - return p; -} - -static NPY_INLINE void PyArray_free_aligned(void *p) -{ - void *base = *(((void **)p) - 1); - PyMem_Free(base); -} - -#endif diff --git a/numpy/random/src/distributions/distributions.c b/numpy/random/src/distributions/distributions.c index 65257ecbf..b382ead0b 100644 --- a/numpy/random/src/distributions/distributions.c +++ b/numpy/random/src/distributions/distributions.c @@ -1,4 +1,4 @@ -#include "distributions.h" +#include "include/distributions.h" #include "ziggurat_constants.h" #include "logfactorial.h" @@ -6,90 +6,52 @@ #include <intrin.h> #endif -/* Random generators for external use */ -float random_float(bitgen_t *bitgen_state) { return next_float(bitgen_state); } - -double random_double(bitgen_t *bitgen_state) { - return next_double(bitgen_state); +/* Inline generators for internal use */ +static NPY_INLINE uint32_t next_uint32(bitgen_t *bitgen_state) { + return bitgen_state->next_uint32(bitgen_state->state); } - -static NPY_INLINE double next_standard_exponential(bitgen_t *bitgen_state) { - return -log(1.0 - next_double(bitgen_state)); +static NPY_INLINE uint64_t next_uint64(bitgen_t *bitgen_state) { + return bitgen_state->next_uint64(bitgen_state->state); } -double random_standard_exponential(bitgen_t *bitgen_state) { - return next_standard_exponential(bitgen_state); +static NPY_INLINE float next_float(bitgen_t *bitgen_state) { + return (next_uint32(bitgen_state) >> 9) * (1.0f / 8388608.0f); } -void random_standard_exponential_fill(bitgen_t *bitgen_state, npy_intp cnt, - double *out) { - npy_intp i; - for (i = 0; i < cnt; i++) { - out[i] = next_standard_exponential(bitgen_state); - } +/* Random generators for external use */ +float random_standard_uniform_f(bitgen_t *bitgen_state) { + return next_float(bitgen_state); } -float random_standard_exponential_f(bitgen_t *bitgen_state) { - return -logf(1.0f - next_float(bitgen_state)); +double random_standard_uniform(bitgen_t *bitgen_state) { + return next_double(bitgen_state); } -void random_double_fill(bitgen_t *bitgen_state, npy_intp cnt, double *out) { +void random_standard_uniform_fill(bitgen_t *bitgen_state, npy_intp cnt, double *out) { npy_intp i; for (i = 0; i < cnt; i++) { out[i] = next_double(bitgen_state); } } -#if 0 -double random_gauss(bitgen_t *bitgen_state) { - if (bitgen_state->has_gauss) { - const double temp = bitgen_state->gauss; - bitgen_state->has_gauss = false; - bitgen_state->gauss = 0.0; - return temp; - } else { - double f, x1, x2, r2; - - do { - x1 = 2.0 * next_double(bitgen_state) - 1.0; - x2 = 2.0 * next_double(bitgen_state) - 1.0; - r2 = x1 * x1 + x2 * x2; - } while (r2 >= 1.0 || r2 == 0.0); - /* Polar method, a more efficient version of the Box-Muller approach. */ - f = sqrt(-2.0 * log(r2) / r2); - /* Keep for next call */ - bitgen_state->gauss = f * x1; - bitgen_state->has_gauss = true; - return f * x2; +void random_standard_uniform_fill_f(bitgen_t *bitgen_state, npy_intp cnt, float *out) { + npy_intp i; + for (i = 0; i < cnt; i++) { + out[i] = next_float(bitgen_state); } } -float random_gauss_f(bitgen_t *bitgen_state) { - if (bitgen_state->has_gauss_f) { - const float temp = bitgen_state->gauss_f; - bitgen_state->has_gauss_f = false; - bitgen_state->gauss_f = 0.0f; - return temp; - } else { - float f, x1, x2, r2; - - do { - x1 = 2.0f * next_float(bitgen_state) - 1.0f; - x2 = 2.0f * next_float(bitgen_state) - 1.0f; - r2 = x1 * x1 + x2 * x2; - } while (r2 >= 1.0 || r2 == 0.0); +double random_standard_exponential(bitgen_t *bitgen_state) { + return -log(1.0 - next_double(bitgen_state)); +} - /* Polar method, a more efficient version of the Box-Muller approach. */ - f = sqrtf(-2.0f * logf(r2) / r2); - /* Keep for next call */ - bitgen_state->gauss_f = f * x1; - bitgen_state->has_gauss_f = true; - return f * x2; +void random_standard_exponential_fill(bitgen_t * bitgen_state, npy_intp cnt, double * out) +{ + npy_intp i; + for (i = 0; i < cnt; i++) { + out[i] = random_standard_exponential(bitgen_state); } } -#endif - -static NPY_INLINE double standard_exponential_zig(bitgen_t *bitgen_state); static double standard_exponential_zig_unlikely(bitgen_t *bitgen_state, uint8_t idx, double x) { @@ -101,11 +63,11 @@ static double standard_exponential_zig_unlikely(bitgen_t *bitgen_state, exp(-x)) { return x; } else { - return standard_exponential_zig(bitgen_state); + return random_standard_exponential_zig(bitgen_state); } } -static NPY_INLINE double standard_exponential_zig(bitgen_t *bitgen_state) { +double random_standard_exponential_zig(bitgen_t *bitgen_state) { uint64_t ri; uint8_t idx; double x; @@ -120,20 +82,26 @@ static NPY_INLINE double standard_exponential_zig(bitgen_t *bitgen_state) { return standard_exponential_zig_unlikely(bitgen_state, idx, x); } -double random_standard_exponential_zig(bitgen_t *bitgen_state) { - return standard_exponential_zig(bitgen_state); +void random_standard_exponential_zig_fill(bitgen_t * bitgen_state, npy_intp cnt, double * out) +{ + npy_intp i; + for (i = 0; i < cnt; i++) { + out[i] = random_standard_exponential_zig(bitgen_state); + } +} + +float random_standard_exponential_f(bitgen_t *bitgen_state) { + return -logf(1.0f - next_float(bitgen_state)); } -void random_standard_exponential_zig_fill(bitgen_t *bitgen_state, npy_intp cnt, - double *out) { +void random_standard_exponential_fill_f(bitgen_t * bitgen_state, npy_intp cnt, float * out) +{ npy_intp i; for (i = 0; i < cnt; i++) { - out[i] = standard_exponential_zig(bitgen_state); + out[i] = random_standard_exponential_f(bitgen_state); } } -static NPY_INLINE float standard_exponential_zig_f(bitgen_t *bitgen_state); - static float standard_exponential_zig_unlikely_f(bitgen_t *bitgen_state, uint8_t idx, float x) { if (idx == 0) { @@ -144,11 +112,11 @@ static float standard_exponential_zig_unlikely_f(bitgen_t *bitgen_state, expf(-x)) { return x; } else { - return standard_exponential_zig_f(bitgen_state); + return random_standard_exponential_zig_f(bitgen_state); } } -static NPY_INLINE float standard_exponential_zig_f(bitgen_t *bitgen_state) { +float random_standard_exponential_zig_f(bitgen_t *bitgen_state) { uint32_t ri; uint8_t idx; float x; @@ -163,11 +131,15 @@ static NPY_INLINE float standard_exponential_zig_f(bitgen_t *bitgen_state) { return standard_exponential_zig_unlikely_f(bitgen_state, idx, x); } -float random_standard_exponential_zig_f(bitgen_t *bitgen_state) { - return standard_exponential_zig_f(bitgen_state); +void random_standard_exponential_zig_fill_f(bitgen_t * bitgen_state, npy_intp cnt, float * out) +{ + npy_intp i; + for (i = 0; i < cnt; i++) { + out[i] = random_standard_exponential_zig_f(bitgen_state); + } } -static NPY_INLINE double next_gauss_zig(bitgen_t *bitgen_state) { +double random_standard_normal(bitgen_t *bitgen_state) { uint64_t r; int sign; uint64_t rabs; @@ -202,18 +174,14 @@ static NPY_INLINE double next_gauss_zig(bitgen_t *bitgen_state) { } } -double random_gauss_zig(bitgen_t *bitgen_state) { - return next_gauss_zig(bitgen_state); -} - -void random_gauss_zig_fill(bitgen_t *bitgen_state, npy_intp cnt, double *out) { +void random_standard_normal_fill(bitgen_t *bitgen_state, npy_intp cnt, double *out) { npy_intp i; for (i = 0; i < cnt; i++) { - out[i] = next_gauss_zig(bitgen_state); + out[i] = random_standard_normal(bitgen_state); } } -float random_gauss_zig_f(bitgen_t *bitgen_state) { +float random_standard_normal_f(bitgen_t *bitgen_state) { uint32_t r; int sign; uint32_t rabs; @@ -247,101 +215,14 @@ float random_gauss_zig_f(bitgen_t *bitgen_state) { } } -/* -static NPY_INLINE double standard_gamma(bitgen_t *bitgen_state, double shape) { - double b, c; - double U, V, X, Y; - - if (shape == 1.0) { - return random_standard_exponential(bitgen_state); - } else if (shape < 1.0) { - for (;;) { - U = next_double(bitgen_state); - V = random_standard_exponential(bitgen_state); - if (U <= 1.0 - shape) { - X = pow(U, 1. / shape); - if (X <= V) { - return X; - } - } else { - Y = -log((1 - U) / shape); - X = pow(1.0 - shape + shape * Y, 1. / shape); - if (X <= (V + Y)) { - return X; - } - } - } - } else { - b = shape - 1. / 3.; - c = 1. / sqrt(9 * b); - for (;;) { - do { - X = random_gauss(bitgen_state); - V = 1.0 + c * X; - } while (V <= 0.0); - - V = V * V * V; - U = next_double(bitgen_state); - if (U < 1.0 - 0.0331 * (X * X) * (X * X)) - return (b * V); - if (log(U) < 0.5 * X * X + b * (1. - V + log(V))) - return (b * V); - } - } -} - -static NPY_INLINE float standard_gamma_float(bitgen_t *bitgen_state, float -shape) { float b, c; float U, V, X, Y; - - if (shape == 1.0f) { - return random_standard_exponential_f(bitgen_state); - } else if (shape < 1.0f) { - for (;;) { - U = next_float(bitgen_state); - V = random_standard_exponential_f(bitgen_state); - if (U <= 1.0f - shape) { - X = powf(U, 1.0f / shape); - if (X <= V) { - return X; - } - } else { - Y = -logf((1.0f - U) / shape); - X = powf(1.0f - shape + shape * Y, 1.0f / shape); - if (X <= (V + Y)) { - return X; - } - } - } - } else { - b = shape - 1.0f / 3.0f; - c = 1.0f / sqrtf(9.0f * b); - for (;;) { - do { - X = random_gauss_f(bitgen_state); - V = 1.0f + c * X; - } while (V <= 0.0f); - - V = V * V * V; - U = next_float(bitgen_state); - if (U < 1.0f - 0.0331f * (X * X) * (X * X)) - return (b * V); - if (logf(U) < 0.5f * X * X + b * (1.0f - V + logf(V))) - return (b * V); - } +void random_standard_normal_fill_f(bitgen_t *bitgen_state, npy_intp cnt, float *out) { + npy_intp i; + for (i = 0; i < cnt; i++) { + out[i] = random_standard_normal_f(bitgen_state); } } - -double random_standard_gamma(bitgen_t *bitgen_state, double shape) { - return standard_gamma(bitgen_state, shape); -} - -float random_standard_gamma_f(bitgen_t *bitgen_state, float shape) { - return standard_gamma_float(bitgen_state, shape); -} -*/ - -static NPY_INLINE double standard_gamma_zig(bitgen_t *bitgen_state, +double random_standard_gamma(bitgen_t *bitgen_state, double shape) { double b, c; double U, V, X, Y; @@ -372,7 +253,7 @@ static NPY_INLINE double standard_gamma_zig(bitgen_t *bitgen_state, c = 1. / sqrt(9 * b); for (;;) { do { - X = random_gauss_zig(bitgen_state); + X = random_standard_normal(bitgen_state); V = 1.0 + c * X; } while (V <= 0.0); @@ -387,7 +268,7 @@ static NPY_INLINE double standard_gamma_zig(bitgen_t *bitgen_state, } } -static NPY_INLINE float standard_gamma_zig_f(bitgen_t *bitgen_state, +float random_standard_gamma_f(bitgen_t *bitgen_state, float shape) { float b, c; float U, V, X, Y; @@ -418,7 +299,7 @@ static NPY_INLINE float standard_gamma_zig_f(bitgen_t *bitgen_state, c = 1.0f / sqrtf(9.0f * b); for (;;) { do { - X = random_gauss_zig_f(bitgen_state); + X = random_standard_normal_f(bitgen_state); V = 1.0f + c * X; } while (V <= 0.0f); @@ -433,14 +314,6 @@ static NPY_INLINE float standard_gamma_zig_f(bitgen_t *bitgen_state, } } -double random_standard_gamma_zig(bitgen_t *bitgen_state, double shape) { - return standard_gamma_zig(bitgen_state, shape); -} - -float random_standard_gamma_zig_f(bitgen_t *bitgen_state, float shape) { - return standard_gamma_zig_f(bitgen_state, shape); -} - int64_t random_positive_int64(bitgen_t *bitgen_state) { return next_uint64(bitgen_state) >> 1; } @@ -470,10 +343,10 @@ uint64_t random_uint(bitgen_t *bitgen_state) { * algorithm comes from SPECFUN by Shanjie Zhang and Jianming Jin and their * book "Computation of Special Functions", 1996, John Wiley & Sons, Inc. * - * If loggam(k+1) is being used to compute log(k!) for an integer k, consider + * If random_loggam(k+1) is being used to compute log(k!) for an integer k, consider * using logfactorial(k) instead. */ -double loggam(double x) { +double random_loggam(double x) { double x0, x2, xp, gl, gl0; RAND_INT_TYPE k, n; @@ -513,12 +386,12 @@ double random_normal(bitgen_t *bitgen_state, double loc, double scale) { } */ -double random_normal_zig(bitgen_t *bitgen_state, double loc, double scale) { - return loc + scale * random_gauss_zig(bitgen_state); +double random_normal(bitgen_t *bitgen_state, double loc, double scale) { + return loc + scale * random_standard_normal(bitgen_state); } double random_exponential(bitgen_t *bitgen_state, double scale) { - return scale * standard_exponential_zig(bitgen_state); + return scale * random_standard_exponential_zig(bitgen_state); } double random_uniform(bitgen_t *bitgen_state, double lower, double range) { @@ -526,11 +399,11 @@ double random_uniform(bitgen_t *bitgen_state, double lower, double range) { } double random_gamma(bitgen_t *bitgen_state, double shape, double scale) { - return scale * random_standard_gamma_zig(bitgen_state, shape); + return scale * random_standard_gamma(bitgen_state, shape); } -float random_gamma_float(bitgen_t *bitgen_state, float shape, float scale) { - return scale * random_standard_gamma_zig_f(bitgen_state, shape); +float random_gamma_f(bitgen_t *bitgen_state, float shape, float scale) { + return scale * random_standard_gamma_f(bitgen_state, shape); } double random_beta(bitgen_t *bitgen_state, double a, double b) { @@ -562,14 +435,14 @@ double random_beta(bitgen_t *bitgen_state, double a, double b) { } } } else { - Ga = random_standard_gamma_zig(bitgen_state, a); - Gb = random_standard_gamma_zig(bitgen_state, b); + Ga = random_standard_gamma(bitgen_state, a); + Gb = random_standard_gamma(bitgen_state, b); return Ga / (Ga + Gb); } } double random_chisquare(bitgen_t *bitgen_state, double df) { - return 2.0 * random_standard_gamma_zig(bitgen_state, df / 2.0); + return 2.0 * random_standard_gamma(bitgen_state, df / 2.0); } double random_f(bitgen_t *bitgen_state, double dfnum, double dfden) { @@ -578,22 +451,22 @@ double random_f(bitgen_t *bitgen_state, double dfnum, double dfden) { } double random_standard_cauchy(bitgen_t *bitgen_state) { - return random_gauss_zig(bitgen_state) / random_gauss_zig(bitgen_state); + return random_standard_normal(bitgen_state) / random_standard_normal(bitgen_state); } double random_pareto(bitgen_t *bitgen_state, double a) { - return exp(standard_exponential_zig(bitgen_state) / a) - 1; + return exp(random_standard_exponential_zig(bitgen_state) / a) - 1; } double random_weibull(bitgen_t *bitgen_state, double a) { if (a == 0.0) { return 0.0; } - return pow(standard_exponential_zig(bitgen_state), 1. / a); + return pow(random_standard_exponential_zig(bitgen_state), 1. / a); } double random_power(bitgen_t *bitgen_state, double a) { - return pow(1 - exp(-standard_exponential_zig(bitgen_state)), 1. / a); + return pow(1 - exp(-random_standard_exponential_zig(bitgen_state)), 1. / a); } double random_laplace(bitgen_t *bitgen_state, double loc, double scale) { @@ -634,7 +507,7 @@ double random_logistic(bitgen_t *bitgen_state, double loc, double scale) { } double random_lognormal(bitgen_t *bitgen_state, double mean, double sigma) { - return exp(random_normal_zig(bitgen_state, mean, sigma)); + return exp(random_normal(bitgen_state, mean, sigma)); } double random_rayleigh(bitgen_t *bitgen_state, double mode) { @@ -644,8 +517,8 @@ double random_rayleigh(bitgen_t *bitgen_state, double mode) { double random_standard_t(bitgen_t *bitgen_state, double df) { double num, denom; - num = random_gauss_zig(bitgen_state); - denom = random_standard_gamma_zig(bitgen_state, df / 2); + num = random_standard_normal(bitgen_state); + denom = random_standard_gamma(bitgen_state, df / 2); return sqrt(df / 2) * num / sqrt(denom); } @@ -699,7 +572,7 @@ static RAND_INT_TYPE random_poisson_ptrs(bitgen_t *bitgen_state, double lam) { /* log(V) == log(0.0) ok here */ /* if U==0.0 so that us==0.0, log is ok since always returns */ if ((log(V) + log(invalpha) - log(a / (us * us) + b)) <= - (-lam + k * loglam - loggam(k + 1))) { + (-lam + k * loglam - random_loggam(k + 1))) { return k; } } @@ -901,8 +774,8 @@ RAND_INT_TYPE random_binomial_inversion(bitgen_t *bitgen_state, RAND_INT_TYPE n, return X; } -RAND_INT_TYPE random_binomial(bitgen_t *bitgen_state, double p, RAND_INT_TYPE n, - binomial_t *binomial) { +int64_t random_binomial(bitgen_t *bitgen_state, double p, int64_t n, + binomial_t *binomial) { double q; if ((n == 0LL) || (p == 0.0f)) @@ -934,7 +807,7 @@ double random_noncentral_chisquare(bitgen_t *bitgen_state, double df, } if (1 < df) { const double Chi2 = random_chisquare(bitgen_state, df - 1); - const double n = random_gauss_zig(bitgen_state) + sqrt(nonc); + const double n = random_standard_normal(bitgen_state) + sqrt(nonc); return Chi2 + n * n; } else { const RAND_INT_TYPE i = random_poisson(bitgen_state, nonc / 2.0); @@ -953,7 +826,7 @@ double random_wald(bitgen_t *bitgen_state, double mean, double scale) { double mu_2l; mu_2l = mean / (2 * scale); - Y = random_gauss_zig(bitgen_state); + Y = random_standard_normal(bitgen_state); Y = mean * Y * Y; X = mean + mu_2l * (Y - sqrt(4 * scale * Y + Y * Y)); U = next_double(bitgen_state); @@ -1092,8 +965,8 @@ RAND_INT_TYPE random_zipf(bitgen_t *bitgen_state, double a) { while (1) { double T, U, V, X; - U = 1.0 - random_double(bitgen_state); - V = random_double(bitgen_state); + U = 1.0 - next_double(bitgen_state); + V = next_double(bitgen_state); X = floor(pow(U, -1.0 / am1)); /* * The real result may be above what can be represented in a signed @@ -1478,7 +1351,7 @@ uint64_t random_bounded_uint64(bitgen_t *bitgen_state, uint64_t off, uint64_t rng, uint64_t mask, bool use_masked) { if (rng == 0) { return off; - } else if (rng < 0xFFFFFFFFUL) { + } else if (rng <= 0xFFFFFFFFUL) { /* Call 32-bit generator if range in 32-bit. */ if (use_masked) { return off + buffered_bounded_masked_uint32(bitgen_state, rng, mask, NULL, @@ -1592,7 +1465,7 @@ void random_bounded_uint64_fill(bitgen_t *bitgen_state, uint64_t off, for (i = 0; i < cnt; i++) { out[i] = off; } - } else if (rng < 0xFFFFFFFFUL) { + } else if (rng <= 0xFFFFFFFFUL) { uint32_t buf = 0; int bcnt = 0; diff --git a/numpy/random/src/distributions/distributions.h b/numpy/random/src/distributions/distributions.h deleted file mode 100644 index c8cdfd20f..000000000 --- a/numpy/random/src/distributions/distributions.h +++ /dev/null @@ -1,205 +0,0 @@ -#ifndef _RANDOMDGEN__DISTRIBUTIONS_H_ -#define _RANDOMDGEN__DISTRIBUTIONS_H_ - -#pragma once -#include <stddef.h> -#include <stdbool.h> -#include <stdint.h> - -#include "Python.h" -#include "numpy/npy_common.h" -#include "numpy/npy_math.h" -#include "numpy/random/bitgen.h" - -/* - * RAND_INT_TYPE is used to share integer generators with RandomState which - * used long in place of int64_t. If changing a distribution that uses - * RAND_INT_TYPE, then the original unmodified copy must be retained for - * use in RandomState by copying to the legacy distributions source file. - */ -#ifdef NP_RANDOM_LEGACY -#define RAND_INT_TYPE long -#define RAND_INT_MAX LONG_MAX -#else -#define RAND_INT_TYPE int64_t -#define RAND_INT_MAX INT64_MAX -#endif - -#ifdef DLL_EXPORT -#define DECLDIR __declspec(dllexport) -#else -#define DECLDIR extern -#endif - -#ifndef MIN -#define MIN(x, y) (((x) < (y)) ? x : y) -#define MAX(x, y) (((x) > (y)) ? x : y) -#endif - -#ifndef M_PI -#define M_PI 3.14159265358979323846264338328 -#endif - -typedef struct s_binomial_t { - int has_binomial; /* !=0: following parameters initialized for binomial */ - double psave; - int64_t nsave; - double r; - double q; - double fm; - int64_t m; - double p1; - double xm; - double xl; - double xr; - double c; - double laml; - double lamr; - double p2; - double p3; - double p4; -} binomial_t; - -/* Inline generators for internal use */ -static NPY_INLINE uint32_t next_uint32(bitgen_t *bitgen_state) { - return bitgen_state->next_uint32(bitgen_state->state); -} - -static NPY_INLINE uint64_t next_uint64(bitgen_t *bitgen_state) { - return bitgen_state->next_uint64(bitgen_state->state); -} - -static NPY_INLINE float next_float(bitgen_t *bitgen_state) { - return (next_uint32(bitgen_state) >> 9) * (1.0f / 8388608.0f); -} - -static NPY_INLINE double next_double(bitgen_t *bitgen_state) { - return bitgen_state->next_double(bitgen_state->state); -} - -DECLDIR double loggam(double x); - -DECLDIR float random_float(bitgen_t *bitgen_state); -DECLDIR double random_double(bitgen_t *bitgen_state); -DECLDIR void random_double_fill(bitgen_t *bitgen_state, npy_intp cnt, double *out); - -DECLDIR int64_t random_positive_int64(bitgen_t *bitgen_state); -DECLDIR int32_t random_positive_int32(bitgen_t *bitgen_state); -DECLDIR int64_t random_positive_int(bitgen_t *bitgen_state); -DECLDIR uint64_t random_uint(bitgen_t *bitgen_state); - -DECLDIR double random_standard_exponential(bitgen_t *bitgen_state); -DECLDIR void random_standard_exponential_fill(bitgen_t *bitgen_state, npy_intp cnt, - double *out); -DECLDIR float random_standard_exponential_f(bitgen_t *bitgen_state); -DECLDIR double random_standard_exponential_zig(bitgen_t *bitgen_state); -DECLDIR void random_standard_exponential_zig_fill(bitgen_t *bitgen_state, - npy_intp cnt, double *out); -DECLDIR float random_standard_exponential_zig_f(bitgen_t *bitgen_state); - -/* -DECLDIR double random_gauss(bitgen_t *bitgen_state); -DECLDIR float random_gauss_f(bitgen_t *bitgen_state); -*/ -DECLDIR double random_gauss_zig(bitgen_t *bitgen_state); -DECLDIR float random_gauss_zig_f(bitgen_t *bitgen_state); -DECLDIR void random_gauss_zig_fill(bitgen_t *bitgen_state, npy_intp cnt, - double *out); - -/* -DECLDIR double random_standard_gamma(bitgen_t *bitgen_state, double shape); -DECLDIR float random_standard_gamma_f(bitgen_t *bitgen_state, float shape); -*/ -DECLDIR double random_standard_gamma_zig(bitgen_t *bitgen_state, double shape); -DECLDIR float random_standard_gamma_zig_f(bitgen_t *bitgen_state, float shape); - -/* -DECLDIR double random_normal(bitgen_t *bitgen_state, double loc, double scale); -*/ -DECLDIR double random_normal_zig(bitgen_t *bitgen_state, double loc, double scale); - -DECLDIR double random_gamma(bitgen_t *bitgen_state, double shape, double scale); -DECLDIR float random_gamma_float(bitgen_t *bitgen_state, float shape, float scale); - -DECLDIR double random_exponential(bitgen_t *bitgen_state, double scale); -DECLDIR double random_uniform(bitgen_t *bitgen_state, double lower, double range); -DECLDIR double random_beta(bitgen_t *bitgen_state, double a, double b); -DECLDIR double random_chisquare(bitgen_t *bitgen_state, double df); -DECLDIR double random_f(bitgen_t *bitgen_state, double dfnum, double dfden); -DECLDIR double random_standard_cauchy(bitgen_t *bitgen_state); -DECLDIR double random_pareto(bitgen_t *bitgen_state, double a); -DECLDIR double random_weibull(bitgen_t *bitgen_state, double a); -DECLDIR double random_power(bitgen_t *bitgen_state, double a); -DECLDIR double random_laplace(bitgen_t *bitgen_state, double loc, double scale); -DECLDIR double random_gumbel(bitgen_t *bitgen_state, double loc, double scale); -DECLDIR double random_logistic(bitgen_t *bitgen_state, double loc, double scale); -DECLDIR double random_lognormal(bitgen_t *bitgen_state, double mean, double sigma); -DECLDIR double random_rayleigh(bitgen_t *bitgen_state, double mode); -DECLDIR double random_standard_t(bitgen_t *bitgen_state, double df); -DECLDIR double random_noncentral_chisquare(bitgen_t *bitgen_state, double df, - double nonc); -DECLDIR double random_noncentral_f(bitgen_t *bitgen_state, double dfnum, - double dfden, double nonc); -DECLDIR double random_wald(bitgen_t *bitgen_state, double mean, double scale); -DECLDIR double random_vonmises(bitgen_t *bitgen_state, double mu, double kappa); -DECLDIR double random_triangular(bitgen_t *bitgen_state, double left, double mode, - double right); - -DECLDIR RAND_INT_TYPE random_poisson(bitgen_t *bitgen_state, double lam); -DECLDIR RAND_INT_TYPE random_negative_binomial(bitgen_t *bitgen_state, double n, - double p); -DECLDIR RAND_INT_TYPE random_binomial(bitgen_t *bitgen_state, double p, RAND_INT_TYPE n, - binomial_t *binomial); -DECLDIR RAND_INT_TYPE random_logseries(bitgen_t *bitgen_state, double p); -DECLDIR RAND_INT_TYPE random_geometric_search(bitgen_t *bitgen_state, double p); -DECLDIR RAND_INT_TYPE random_geometric_inversion(bitgen_t *bitgen_state, double p); -DECLDIR RAND_INT_TYPE random_geometric(bitgen_t *bitgen_state, double p); -DECLDIR RAND_INT_TYPE random_zipf(bitgen_t *bitgen_state, double a); -DECLDIR int64_t random_hypergeometric(bitgen_t *bitgen_state, - int64_t good, int64_t bad, int64_t sample); - -DECLDIR uint64_t random_interval(bitgen_t *bitgen_state, uint64_t max); - -/* Generate random uint64 numbers in closed interval [off, off + rng]. */ -DECLDIR uint64_t random_bounded_uint64(bitgen_t *bitgen_state, uint64_t off, - uint64_t rng, uint64_t mask, - bool use_masked); - -/* Generate random uint32 numbers in closed interval [off, off + rng]. */ -DECLDIR uint32_t random_buffered_bounded_uint32(bitgen_t *bitgen_state, - uint32_t off, uint32_t rng, - uint32_t mask, bool use_masked, - int *bcnt, uint32_t *buf); -DECLDIR uint16_t random_buffered_bounded_uint16(bitgen_t *bitgen_state, - uint16_t off, uint16_t rng, - uint16_t mask, bool use_masked, - int *bcnt, uint32_t *buf); -DECLDIR uint8_t random_buffered_bounded_uint8(bitgen_t *bitgen_state, uint8_t off, - uint8_t rng, uint8_t mask, - bool use_masked, int *bcnt, - uint32_t *buf); -DECLDIR npy_bool random_buffered_bounded_bool(bitgen_t *bitgen_state, npy_bool off, - npy_bool rng, npy_bool mask, - bool use_masked, int *bcnt, - uint32_t *buf); - -DECLDIR void random_bounded_uint64_fill(bitgen_t *bitgen_state, uint64_t off, - uint64_t rng, npy_intp cnt, - bool use_masked, uint64_t *out); -DECLDIR void random_bounded_uint32_fill(bitgen_t *bitgen_state, uint32_t off, - uint32_t rng, npy_intp cnt, - bool use_masked, uint32_t *out); -DECLDIR void random_bounded_uint16_fill(bitgen_t *bitgen_state, uint16_t off, - uint16_t rng, npy_intp cnt, - bool use_masked, uint16_t *out); -DECLDIR void random_bounded_uint8_fill(bitgen_t *bitgen_state, uint8_t off, - uint8_t rng, npy_intp cnt, - bool use_masked, uint8_t *out); -DECLDIR void random_bounded_bool_fill(bitgen_t *bitgen_state, npy_bool off, - npy_bool rng, npy_intp cnt, - bool use_masked, npy_bool *out); - -DECLDIR void random_multinomial(bitgen_t *bitgen_state, RAND_INT_TYPE n, RAND_INT_TYPE *mnix, - double *pix, npy_intp d, binomial_t *binomial); - -#endif diff --git a/numpy/random/src/distributions/random_hypergeometric.c b/numpy/random/src/distributions/random_hypergeometric.c index 59a3a4b9b..da5ea9c68 100644 --- a/numpy/random/src/distributions/random_hypergeometric.c +++ b/numpy/random/src/distributions/random_hypergeometric.c @@ -1,6 +1,6 @@ -#include <stdint.h> -#include "distributions.h" +#include "include/distributions.h" #include "logfactorial.h" +#include <stdint.h> /* * Generate a sample from the hypergeometric distribution. @@ -188,8 +188,8 @@ static int64_t hypergeometric_hrua(bitgen_t *bitgen_state, while (1) { double U, V, X, T; double gp; - U = random_double(bitgen_state); - V = random_double(bitgen_state); // "U star" in Stadlober (1989) + U = next_double(bitgen_state); + V = next_double(bitgen_state); // "U star" in Stadlober (1989) X = a + h*(V - 0.5) / U; // fast rejection: diff --git a/numpy/random/src/distributions/random_mvhg_count.c b/numpy/random/src/distributions/random_mvhg_count.c new file mode 100644 index 000000000..9c0cc045d --- /dev/null +++ b/numpy/random/src/distributions/random_mvhg_count.c @@ -0,0 +1,131 @@ +#include <stdint.h> +#include <stdlib.h> +#include <stdbool.h> + +#include "include/distributions.h" + +/* + * random_mvhg_count + * + * Draw variates from the multivariate hypergeometric distribution-- + * the "count" algorithm. + * + * Parameters + * ---------- + * bitgen_t *bitgen_state + * Pointer to a `bitgen_t` instance. + * int64_t total + * The sum of the values in the array `colors`. (This is redundant + * information, but we know the caller has already computed it, so + * we might as well use it.) + * size_t num_colors + * The length of the `colors` array. + * int64_t *colors + * The array of colors (i.e. the number of each type in the collection + * from which the random variate is drawn). + * int64_t nsample + * The number of objects drawn without replacement for each variate. + * `nsample` must not exceed sum(colors). This condition is not checked; + * it is assumed that the caller has already validated the value. + * size_t num_variates + * The number of variates to be produced and put in the array + * pointed to by `variates`. One variate is a vector of length + * `num_colors`, so the array pointed to by `variates` must have length + * `num_variates * num_colors`. + * int64_t *variates + * The array that will hold the result. It must have length + * `num_variates * num_colors`. + * The array is not initialized in the function; it is expected that the + * array has been initialized with zeros when the function is called. + * + * Notes + * ----- + * The "count" algorithm for drawing one variate is roughly equivalent to the + * following numpy code: + * + * choices = np.repeat(np.arange(len(colors)), colors) + * selection = np.random.choice(choices, nsample, replace=False) + * variate = np.bincount(selection, minlength=len(colors)) + * + * This function uses a temporary array with length sum(colors). + * + * Assumptions on the arguments (not checked in the function): + * * colors[k] >= 0 for k in range(num_colors) + * * total = sum(colors) + * * 0 <= nsample <= total + * * the product total * sizeof(size_t) does not exceed SIZE_MAX + * * the product num_variates * num_colors does not overflow + */ + +int random_mvhg_count(bitgen_t *bitgen_state, + int64_t total, + size_t num_colors, int64_t *colors, + int64_t nsample, + size_t num_variates, int64_t *variates) +{ + size_t *choices; + bool more_than_half; + + if ((total == 0) || (nsample == 0) || (num_variates == 0)) { + // Nothing to do. + return 0; + } + + choices = malloc(total * (sizeof *choices)); + if (choices == NULL) { + return -1; + } + + /* + * If colors contains, for example, [3 2 5], then choices + * will contain [0 0 0 1 1 2 2 2 2 2]. + */ + for (size_t i = 0, k = 0; i < num_colors; ++i) { + for (int64_t j = 0; j < colors[i]; ++j) { + choices[k] = i; + ++k; + } + } + + more_than_half = nsample > (total / 2); + if (more_than_half) { + nsample = total - nsample; + } + + for (size_t i = 0; i < num_variates * num_colors; i += num_colors) { + /* + * Fisher-Yates shuffle, but only loop through the first + * `nsample` entries of `choices`. After the loop, + * choices[:nsample] contains a random sample from the + * the full array. + */ + for (size_t j = 0; j < (size_t) nsample; ++j) { + size_t tmp, k; + // Note: nsample is not greater than total, so there is no danger + // of integer underflow in `(size_t) total - j - 1`. + k = j + (size_t) random_interval(bitgen_state, + (size_t) total - j - 1); + tmp = choices[k]; + choices[k] = choices[j]; + choices[j] = tmp; + } + /* + * Count the number of occurrences of each value in choices[:nsample]. + * The result, stored in sample[i:i+num_colors], is the sample from + * the multivariate hypergeometric distribution. + */ + for (size_t j = 0; j < (size_t) nsample; ++j) { + variates[i + choices[j]] += 1; + } + + if (more_than_half) { + for (size_t k = 0; k < num_colors; ++k) { + variates[i + k] = colors[k] - variates[i + k]; + } + } + } + + free(choices); + + return 0; +} diff --git a/numpy/random/src/distributions/random_mvhg_marginals.c b/numpy/random/src/distributions/random_mvhg_marginals.c new file mode 100644 index 000000000..301a4acad --- /dev/null +++ b/numpy/random/src/distributions/random_mvhg_marginals.c @@ -0,0 +1,138 @@ +#include <stdint.h> +#include <stddef.h> +#include <stdbool.h> +#include <math.h> + +#include "include/distributions.h" +#include "logfactorial.h" + + +/* + * random_mvhg_marginals + * + * Draw samples from the multivariate hypergeometric distribution-- + * the "marginals" algorithm. + * + * This version generates the sample by iteratively calling + * hypergeometric() (the univariate hypergeometric distribution). + * + * Parameters + * ---------- + * bitgen_t *bitgen_state + * Pointer to a `bitgen_t` instance. + * int64_t total + * The sum of the values in the array `colors`. (This is redundant + * information, but we know the caller has already computed it, so + * we might as well use it.) + * size_t num_colors + * The length of the `colors` array. The functions assumes + * num_colors > 0. + * int64_t *colors + * The array of colors (i.e. the number of each type in the collection + * from which the random variate is drawn). + * int64_t nsample + * The number of objects drawn without replacement for each variate. + * `nsample` must not exceed sum(colors). This condition is not checked; + * it is assumed that the caller has already validated the value. + * size_t num_variates + * The number of variates to be produced and put in the array + * pointed to by `variates`. One variate is a vector of length + * `num_colors`, so the array pointed to by `variates` must have length + * `num_variates * num_colors`. + * int64_t *variates + * The array that will hold the result. It must have length + * `num_variates * num_colors`. + * The array is not initialized in the function; it is expected that the + * array has been initialized with zeros when the function is called. + * + * Notes + * ----- + * Here's an example that demonstrates the idea of this algorithm. + * + * Suppose the urn contains red, green, blue and yellow marbles. + * Let nred be the number of red marbles, and define the quantities for + * the other colors similarly. The total number of marbles is + * + * total = nred + ngreen + nblue + nyellow. + * + * To generate a sample using rk_hypergeometric: + * + * red_sample = hypergeometric(ngood=nred, nbad=total - nred, + * nsample=nsample) + * + * This gives us the number of red marbles in the sample. The number of + * marbles in the sample that are *not* red is nsample - red_sample. + * To figure out the distribution of those marbles, we again use + * rk_hypergeometric: + * + * green_sample = hypergeometric(ngood=ngreen, + * nbad=total - nred - ngreen, + * nsample=nsample - red_sample) + * + * Similarly, + * + * blue_sample = hypergeometric( + * ngood=nblue, + * nbad=total - nred - ngreen - nblue, + * nsample=nsample - red_sample - green_sample) + * + * Finally, + * + * yellow_sample = total - (red_sample + green_sample + blue_sample). + * + * The above sequence of steps is implemented as a loop for an arbitrary + * number of colors in the innermost loop in the code below. `remaining` + * is the value passed to `nbad`; it is `total - colors[0]` in the first + * call to random_hypergeometric(), and then decreases by `colors[j]` in + * each iteration. `num_to_sample` is the `nsample` argument. It + * starts at this function's `nsample` input, and is decreased by the + * result of the call to random_hypergeometric() in each iteration. + * + * Assumptions on the arguments (not checked in the function): + * * colors[k] >= 0 for k in range(num_colors) + * * total = sum(colors) + * * 0 <= nsample <= total + * * the product num_variates * num_colors does not overflow + */ + +void random_mvhg_marginals(bitgen_t *bitgen_state, + int64_t total, + size_t num_colors, int64_t *colors, + int64_t nsample, + size_t num_variates, int64_t *variates) +{ + bool more_than_half; + + if ((total == 0) || (nsample == 0) || (num_variates == 0)) { + // Nothing to do. + return; + } + + more_than_half = nsample > (total / 2); + if (more_than_half) { + nsample = total - nsample; + } + + for (size_t i = 0; i < num_variates * num_colors; i += num_colors) { + int64_t num_to_sample = nsample; + int64_t remaining = total; + for (size_t j = 0; (num_to_sample > 0) && (j + 1 < num_colors); ++j) { + int64_t r; + remaining -= colors[j]; + r = random_hypergeometric(bitgen_state, + colors[j], remaining, num_to_sample); + variates[i + j] = r; + num_to_sample -= r; + } + + if (num_to_sample > 0) { + variates[i + num_colors - 1] = num_to_sample; + } + + if (more_than_half) { + for (size_t k = 0; k < num_colors; ++k) { + variates[i + k] = colors[k] - variates[i + k]; + } + } + } +} diff --git a/numpy/random/src/entropy/entropy.c b/numpy/random/src/entropy/entropy.c deleted file mode 100644 index eaca37a9c..000000000 --- a/numpy/random/src/entropy/entropy.c +++ /dev/null @@ -1,114 +0,0 @@ -#include <stddef.h> -#include <stdio.h> -#include <stdlib.h> -#include <string.h> - -#include "entropy.h" -#ifdef _WIN32 -/* Windows */ -#include <sys/timeb.h> -#include <time.h> -#include <windows.h> - -#include <wincrypt.h> -#else -/* Unix */ -#include <sys/time.h> -#include <time.h> -#include <unistd.h> -#include <fcntl.h> -#endif - -bool entropy_getbytes(void *dest, size_t size) { -#ifndef _WIN32 - - int fd = open("/dev/urandom", O_RDONLY); - if (fd < 0) - return false; - ssize_t sz = read(fd, dest, size); - if ((sz < 0) || ((size_t)sz < size)) - return false; - return close(fd) == 0; - -#else - - HCRYPTPROV hCryptProv; - BOOL done; - - if (!CryptAcquireContext(&hCryptProv, NULL, NULL, PROV_RSA_FULL, - CRYPT_VERIFYCONTEXT) || - !hCryptProv) { - return true; - } - done = CryptGenRandom(hCryptProv, (DWORD)size, (unsigned char *)dest); - CryptReleaseContext(hCryptProv, 0); - if (!done) { - return false; - } - - return true; -#endif -} - -/* Thomas Wang 32/64 bits integer hash function */ -uint32_t entropy_hash_32(uint32_t key) { - key += ~(key << 15); - key ^= (key >> 10); - key += (key << 3); - key ^= (key >> 6); - key += ~(key << 11); - key ^= (key >> 16); - return key; -} - -uint64_t entropy_hash_64(uint64_t key) { - key = (~key) + (key << 21); // key = (key << 21) - key - 1; - key = key ^ (key >> 24); - key = (key + (key << 3)) + (key << 8); // key * 265 - key = key ^ (key >> 14); - key = (key + (key << 2)) + (key << 4); // key * 21 - key = key ^ (key >> 28); - key = key + (key << 31); - return key; -} - -uint32_t entropy_randombytes(void) { - -#ifndef _WIN32 - struct timeval tv; - gettimeofday(&tv, NULL); - return entropy_hash_32(getpid()) ^ entropy_hash_32(tv.tv_sec) ^ - entropy_hash_32(tv.tv_usec) ^ entropy_hash_32(clock()); -#else - uint32_t out = 0; - int64_t counter; - struct _timeb tv; - _ftime_s(&tv); - out = entropy_hash_32(GetCurrentProcessId()) ^ - entropy_hash_32((uint32_t)tv.time) ^ entropy_hash_32(tv.millitm) ^ - entropy_hash_32(clock()); - if (QueryPerformanceCounter((LARGE_INTEGER *)&counter) != 0) - out ^= entropy_hash_32((uint32_t)(counter & 0xffffffff)); - return out; -#endif -} - -bool entropy_fallback_getbytes(void *dest, size_t size) { - int hashes = (int)size; - uint32_t *hash = malloc(hashes * sizeof(uint32_t)); - int i; - for (i = 0; i < hashes; i++) { - hash[i] = entropy_randombytes(); - } - memcpy(dest, (void *)hash, size); - free(hash); - return true; -} - -void entropy_fill(void *dest, size_t size) { - bool success; - success = entropy_getbytes(dest, size); - if (!success) { - entropy_fallback_getbytes(dest, size); - } -} diff --git a/numpy/random/src/entropy/entropy.h b/numpy/random/src/entropy/entropy.h deleted file mode 100644 index f00caf61d..000000000 --- a/numpy/random/src/entropy/entropy.h +++ /dev/null @@ -1,14 +0,0 @@ -#ifndef _RANDOMDGEN__ENTROPY_H_ -#define _RANDOMDGEN__ENTROPY_H_ - -#include <stddef.h> -#include <stdbool.h> -#include <stdint.h> - -extern void entropy_fill(void *dest, size_t size); - -extern bool entropy_getbytes(void *dest, size_t size); - -extern bool entropy_fallback_getbytes(void *dest, size_t size); - -#endif diff --git a/numpy/random/src/legacy/legacy-distributions.c b/numpy/random/src/legacy/legacy-distributions.c index 4741a0352..fd067fe8d 100644 --- a/numpy/random/src/legacy/legacy-distributions.c +++ b/numpy/random/src/legacy/legacy-distributions.c @@ -1,4 +1,4 @@ -#include "legacy-distributions.h" +#include "include/legacy-distributions.h" static NPY_INLINE double legacy_double(aug_bitgen_t *aug_state) { @@ -215,6 +215,37 @@ double legacy_exponential(aug_bitgen_t *aug_state, double scale) { } +static RAND_INT_TYPE legacy_random_binomial_original(bitgen_t *bitgen_state, + double p, + RAND_INT_TYPE n, + binomial_t *binomial) { + double q; + + if (p <= 0.5) { + if (p * n <= 30.0) { + return random_binomial_inversion(bitgen_state, n, p, binomial); + } else { + return random_binomial_btpe(bitgen_state, n, p, binomial); + } + } else { + q = 1.0 - p; + if (q * n <= 30.0) { + return n - random_binomial_inversion(bitgen_state, n, q, binomial); + } else { + return n - random_binomial_btpe(bitgen_state, n, q, binomial); + } + } +} + + +int64_t legacy_random_binomial(bitgen_t *bitgen_state, double p, + int64_t n, binomial_t *binomial) { + return (int64_t) legacy_random_binomial_original(bitgen_state, p, + (RAND_INT_TYPE) n, + binomial); +} + + static RAND_INT_TYPE random_hypergeometric_hyp(bitgen_t *bitgen_state, RAND_INT_TYPE good, RAND_INT_TYPE bad, @@ -263,8 +294,8 @@ static RAND_INT_TYPE random_hypergeometric_hrua(bitgen_t *bitgen_state, d7 = sqrt((double)(popsize - m) * sample * d4 * d5 / (popsize - 1) + 0.5); d8 = D1 * d7 + D2; d9 = (RAND_INT_TYPE)floor((double)(m + 1) * (mingoodbad + 1) / (popsize + 2)); - d10 = (loggam(d9 + 1) + loggam(mingoodbad - d9 + 1) + loggam(m - d9 + 1) + - loggam(maxgoodbad - m + d9 + 1)); + d10 = (random_loggam(d9 + 1) + random_loggam(mingoodbad - d9 + 1) + + random_loggam(m - d9 + 1) + random_loggam(maxgoodbad - m + d9 + 1)); d11 = MIN(MIN(m, mingoodbad) + 1.0, floor(d6 + 16 * d7)); /* 16 for 16-decimal-digit precision in D1 and D2 */ @@ -278,8 +309,8 @@ static RAND_INT_TYPE random_hypergeometric_hrua(bitgen_t *bitgen_state, continue; Z = (RAND_INT_TYPE)floor(W); - T = d10 - (loggam(Z + 1) + loggam(mingoodbad - Z + 1) + loggam(m - Z + 1) + - loggam(maxgoodbad - m + Z + 1)); + T = d10 - (random_loggam(Z + 1) + random_loggam(mingoodbad - Z + 1) + + random_loggam(m - Z + 1) + random_loggam(maxgoodbad - m + Z + 1)); /* fast acceptance: */ if ((X * (4.0 - X) - 3.0) <= T) diff --git a/numpy/random/src/legacy/legacy-distributions.h b/numpy/random/src/legacy/legacy-distributions.h deleted file mode 100644 index 005c4e5d2..000000000 --- a/numpy/random/src/legacy/legacy-distributions.h +++ /dev/null @@ -1,52 +0,0 @@ -#ifndef _RANDOMDGEN__DISTRIBUTIONS_LEGACY_H_ -#define _RANDOMDGEN__DISTRIBUTIONS_LEGACY_H_ - - -#include "../distributions/distributions.h" - -typedef struct aug_bitgen { - bitgen_t *bit_generator; - int has_gauss; - double gauss; -} aug_bitgen_t; - -extern double legacy_gauss(aug_bitgen_t *aug_state); -extern double legacy_standard_exponential(aug_bitgen_t *aug_state); -extern double legacy_pareto(aug_bitgen_t *aug_state, double a); -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_pareto(aug_bitgen_t *aug_state, double a); -extern double legacy_weibull(aug_bitgen_t *aug_state, double a); -extern double legacy_chisquare(aug_bitgen_t *aug_state, double df); -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, - double dfden, double nonc); -extern double legacy_wald(aug_bitgen_t *aug_state, double mean, double scale); -extern double legacy_lognormal(aug_bitgen_t *aug_state, double mean, - double sigma); -extern double legacy_standard_t(aug_bitgen_t *aug_state, double df); -extern int64_t legacy_negative_binomial(aug_bitgen_t *aug_state, double n, - double p); -extern double legacy_standard_cauchy(aug_bitgen_t *state); -extern double legacy_beta(aug_bitgen_t *aug_state, double a, double b); -extern double legacy_f(aug_bitgen_t *aug_state, double dfnum, double dfden); -extern double legacy_normal(aug_bitgen_t *aug_state, double loc, double scale); -extern double legacy_standard_gamma(aug_bitgen_t *aug_state, double shape); -extern double legacy_exponential(aug_bitgen_t *aug_state, double scale); -extern int64_t legacy_negative_binomial(aug_bitgen_t *aug_state, double n, - double p); -extern int64_t legacy_random_hypergeometric(bitgen_t *bitgen_state, - int64_t good, int64_t bad, - int64_t sample); -extern int64_t legacy_random_logseries(bitgen_t *bitgen_state, double p); -extern int64_t legacy_random_poisson(bitgen_t *bitgen_state, double lam); -extern int64_t legacy_random_zipf(bitgen_t *bitgen_state, double a); -extern int64_t legacy_random_geometric(bitgen_t *bitgen_state, double p); -void legacy_random_multinomial(bitgen_t *bitgen_state, RAND_INT_TYPE n, - RAND_INT_TYPE *mnix, double *pix, npy_intp d, - binomial_t *binomial); - -#endif diff --git a/numpy/random/src/philox/philox.h b/numpy/random/src/philox/philox.h index 309d89eae..c72424a97 100644 --- a/numpy/random/src/philox/philox.h +++ b/numpy/random/src/philox/philox.h @@ -1,8 +1,8 @@ #ifndef _RANDOMDGEN__PHILOX_H_ #define _RANDOMDGEN__PHILOX_H_ -#include <inttypes.h> #include "numpy/npy_common.h" +#include <inttypes.h> #define PHILOX_BUFFER_SIZE 4L diff --git a/numpy/random/src/sfc64/sfc64.h b/numpy/random/src/sfc64/sfc64.h index 6674ae69c..75c4118d3 100644 --- a/numpy/random/src/sfc64/sfc64.h +++ b/numpy/random/src/sfc64/sfc64.h @@ -1,11 +1,11 @@ #ifndef _RANDOMDGEN__SFC64_H_ #define _RANDOMDGEN__SFC64_H_ +#include "numpy/npy_common.h" #include <inttypes.h> #ifdef _WIN32 #include <stdlib.h> #endif -#include "numpy/npy_common.h" typedef struct s_sfc64_state { uint64_t s[4]; |