diff options
| -rw-r--r-- | numpy/random/mtrand/distributions.c | 21 |
1 files changed, 14 insertions, 7 deletions
diff --git a/numpy/random/mtrand/distributions.c b/numpy/random/mtrand/distributions.c index a7acbd91b..0e64913b7 100644 --- a/numpy/random/mtrand/distributions.c +++ b/numpy/random/mtrand/distributions.c @@ -676,16 +676,23 @@ long rk_zipf(rk_state *state, double a) { double T, U, V; long X; - double b; - - b = pow(2.0, a-1.0); + double am1, b; + + am1 = a - 1.0; + b = pow(2.0, am1); do { - U = rk_double(state); + U = 1.0-rk_double(state); V = rk_double(state); - X = (long)floor(pow(U, -1.0/(a-1.0))); - T = pow(1.0 + 1.0/X, a-1.0); - } while ((V *X*(T-1.0)/(b-1.0)) > (T/b)); + X = (long)floor(pow(U, -1.0/am1)); + /* The real result may be above what can be represented in a signed + * long. It will get casted to -sys.maxint-1. Since this is + * a straightforward rejection algorithm, we can just reject this value + * in the rejection condition below. This function then models a Zipf + * distribution truncated to sys.maxint. + */ + T = pow(1.0 + 1.0/X, am1); + } while (((V*X*(T-1.0)/(b-1.0)) > (T/b)) || X < 1); return X; } |