diff options
Diffstat (limited to 'numpy/random')
-rw-r--r-- | numpy/random/include/distributions.h | 49 | ||||
-rw-r--r-- | numpy/random/src/distributions/distributions.c | 18 | ||||
-rw-r--r-- | numpy/random/src/legacy/legacy-distributions.c | 8 |
3 files changed, 36 insertions, 39 deletions
diff --git a/numpy/random/include/distributions.h b/numpy/random/include/distributions.h index b5439ff33..fb69c7d2c 100644 --- a/numpy/random/include/distributions.h +++ b/numpy/random/include/distributions.h @@ -59,25 +59,6 @@ typedef struct s_binomial_t { 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_standard_uniform_f(bitgen_t *bitgen_state); DECLDIR double random_standard_uniform(bitgen_t *bitgen_state); DECLDIR void random_standard_uniform_fill(bitgen_t *, npy_intp, double *); @@ -135,27 +116,16 @@ DECLDIR double random_triangular(bitgen_t *bitgen_state, double left, double mod 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_btpe(bitgen_t *bitgen_state, - RAND_INT_TYPE n, - double p, - binomial_t *binomial); -DECLDIR RAND_INT_TYPE random_binomial_inversion(bitgen_t *bitgen_state, - RAND_INT_TYPE n, - double p, - binomial_t *binomial); + double p); + DECLDIR int64_t random_binomial(bitgen_t *bitgen_state, double p, int64_t 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]. */ @@ -200,4 +170,19 @@ DECLDIR void random_bounded_bool_fill(bitgen_t *bitgen_state, npy_bool off, DECLDIR void random_multinomial(bitgen_t *bitgen_state, RAND_INT_TYPE n, RAND_INT_TYPE *mnix, double *pix, npy_intp d, binomial_t *binomial); +/* Common to legacy-distributions.c and distributions.c but not exported */ + +RAND_INT_TYPE random_binomial_btpe(bitgen_t *bitgen_state, + RAND_INT_TYPE n, + double p, + binomial_t *binomial); +RAND_INT_TYPE random_binomial_inversion(bitgen_t *bitgen_state, + RAND_INT_TYPE n, + double p, + binomial_t *binomial); +double random_loggam(double x); +static NPY_INLINE double next_double(bitgen_t *bitgen_state) { + return bitgen_state->next_double(bitgen_state->state); +} + #endif diff --git a/numpy/random/src/distributions/distributions.c b/numpy/random/src/distributions/distributions.c index d1475b685..b382ead0b 100644 --- a/numpy/random/src/distributions/distributions.c +++ b/numpy/random/src/distributions/distributions.c @@ -6,6 +6,18 @@ #include <intrin.h> #endif +/* 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); +} + /* Random generators for external use */ float random_standard_uniform_f(bitgen_t *bitgen_state) { return next_float(bitgen_state); @@ -331,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; @@ -560,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; } } diff --git a/numpy/random/src/legacy/legacy-distributions.c b/numpy/random/src/legacy/legacy-distributions.c index 404c3e301..fd067fe8d 100644 --- a/numpy/random/src/legacy/legacy-distributions.c +++ b/numpy/random/src/legacy/legacy-distributions.c @@ -294,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 */ @@ -309,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) |