From cced5101daf6624b40999637bbb0996f11e3a691 Mon Sep 17 00:00:00 2001 From: Robert Kern Date: Wed, 9 Jul 2008 01:55:33 +0000 Subject: BUG: Make sure the Zipf results are within the allowable range. --- numpy/random/mtrand/distributions.c | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) (limited to 'numpy') 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; } -- cgit v1.2.1