diff options
Diffstat (limited to 'numpy/random/src')
| -rw-r--r-- | numpy/random/src/distributions/distributions.c | 22 |
1 files changed, 11 insertions, 11 deletions
diff --git a/numpy/random/src/distributions/distributions.c b/numpy/random/src/distributions/distributions.c index f47c54a53..3df819d9d 100644 --- a/numpy/random/src/distributions/distributions.c +++ b/numpy/random/src/distributions/distributions.c @@ -47,7 +47,7 @@ static double standard_exponential_unlikely(bitgen_t *bitgen_state, uint8_t idx, double x) { if (idx == 0) { /* Switch to 1.0 - U to avoid log(0.0), see GH 13361 */ - return ziggurat_exp_r - log(1.0 - next_double(bitgen_state)); + return ziggurat_exp_r - npy_log1p(-next_double(bitgen_state)); } else if ((fe_double[idx - 1] - fe_double[idx]) * next_double(bitgen_state) + fe_double[idx] < exp(-x)) { @@ -84,7 +84,7 @@ static float standard_exponential_unlikely_f(bitgen_t *bitgen_state, uint8_t idx, float x) { if (idx == 0) { /* Switch to 1.0 - U to avoid log(0.0), see GH 13361 */ - return ziggurat_exp_r_f - logf(1.0f - next_float(bitgen_state)); + return ziggurat_exp_r_f - npy_log1pf(-next_float(bitgen_state)); } else if ((fe_float[idx - 1] - fe_float[idx]) * next_float(bitgen_state) + fe_float[idx] < expf(-x)) { @@ -121,7 +121,7 @@ void random_standard_exponential_inv_fill(bitgen_t * bitgen_state, npy_intp cnt, { npy_intp i; for (i = 0; i < cnt; i++) { - out[i] = -log(1.0 - next_double(bitgen_state)); + out[i] = -npy_log1p(-next_double(bitgen_state)); } } @@ -129,7 +129,7 @@ void random_standard_exponential_inv_fill_f(bitgen_t * bitgen_state, npy_intp cn { npy_intp i; for (i = 0; i < cnt; i++) { - out[i] = -log(1.0 - next_float(bitgen_state)); + out[i] = -npy_log1p(-next_float(bitgen_state)); } } @@ -155,8 +155,8 @@ double random_standard_normal(bitgen_t *bitgen_state) { if (idx == 0) { for (;;) { /* Switch to 1.0 - U to avoid log(0.0), see GH 13361 */ - xx = -ziggurat_nor_inv_r * log(1.0 - next_double(bitgen_state)); - yy = -log(1.0 - next_double(bitgen_state)); + xx = -ziggurat_nor_inv_r * npy_log1p(-next_double(bitgen_state)); + yy = -npy_log1p(-next_double(bitgen_state)); if (yy + yy > xx * xx) return ((rabs >> 8) & 0x1) ? -(ziggurat_nor_r + xx) : ziggurat_nor_r + xx; @@ -196,8 +196,8 @@ float random_standard_normal_f(bitgen_t *bitgen_state) { if (idx == 0) { for (;;) { /* Switch to 1.0 - U to avoid log(0.0), see GH 13361 */ - xx = -ziggurat_nor_inv_r_f * logf(1.0f - next_float(bitgen_state)); - yy = -logf(1.0f - next_float(bitgen_state)); + xx = -ziggurat_nor_inv_r_f * npy_log1pf(-next_float(bitgen_state)); + yy = -npy_log1pf(-next_float(bitgen_state)); if (yy + yy > xx * xx) return ((rabs >> 8) & 0x1) ? -(ziggurat_nor_r_f + xx) : ziggurat_nor_r_f + xx; @@ -508,7 +508,7 @@ double random_lognormal(bitgen_t *bitgen_state, double mean, double sigma) { } double random_rayleigh(bitgen_t *bitgen_state, double mode) { - return mode * sqrt(-2.0 * log(1.0 - next_double(bitgen_state))); + return mode * sqrt(-2.0 * npy_log1p(-next_double(bitgen_state))); } double random_standard_t(bitgen_t *bitgen_state, double df) { @@ -916,7 +916,7 @@ RAND_INT_TYPE random_logseries(bitgen_t *bitgen_state, double p) { double q, r, U, V; RAND_INT_TYPE result; - r = log(1.0 - p); + r = npy_log1p(-p); while (1) { V = next_double(bitgen_state); @@ -958,7 +958,7 @@ RAND_INT_TYPE random_geometric_search(bitgen_t *bitgen_state, double p) { } RAND_INT_TYPE random_geometric_inversion(bitgen_t *bitgen_state, double p) { - return (RAND_INT_TYPE)ceil(log(1.0 - next_double(bitgen_state)) / log(1.0 - p)); + return (RAND_INT_TYPE)ceil(npy_log1p(-next_double(bitgen_state)) / npy_log1p(-p)); } RAND_INT_TYPE random_geometric(bitgen_t *bitgen_state, double p) { |