diff options
author | Charles Harris <charlesr.harris@gmail.com> | 2017-10-04 12:22:14 -0600 |
---|---|---|
committer | Charles Harris <charlesr.harris@gmail.com> | 2017-10-05 09:38:07 -0600 |
commit | e6932a8266ca2e1d1fdc7ec78e9254b6563b0e39 (patch) | |
tree | e5155de17ef56fc5005073286ceaa1c2db5aeca5 | |
parent | 8989a31d6ddaaa1f7df0e5da9bdee88cd4900a66 (diff) | |
download | python-numpy-e6932a8266ca2e1d1fdc7ec78e9254b6563b0e39.tar.gz python-numpy-e6932a8266ca2e1d1fdc7ec78e9254b6563b0e39.tar.bz2 python-numpy-e6932a8266ca2e1d1fdc7ec78e9254b6563b0e39.zip |
MAINT: Refactor rk_zipf to be more efficient.
-rw-r--r-- | numpy/random/mtrand/distributions.c | 24 |
1 files changed, 12 insertions, 12 deletions
diff --git a/numpy/random/mtrand/distributions.c b/numpy/random/mtrand/distributions.c index a455140f8..b7e157915 100644 --- a/numpy/random/mtrand/distributions.c +++ b/numpy/random/mtrand/distributions.c @@ -720,31 +720,31 @@ double rk_wald(rk_state *state, double mean, double scale) long rk_zipf(rk_state *state, double a) { - double T, U, V, X; double am1, b; am1 = a - 1.0; b = pow(2.0, am1); - T = 0.0; - do - { - U = 1.0-rk_double(state); + while (1) { + double T, U, V, X; + + U = 1.0 - rk_double(state); V = rk_double(state); X = floor(pow(U, -1.0/am1)); - /* The real result may be above what can be represented in a signed + /* + * The real result may be above what can be represented in a signed * long. Since this is a straightforward rejection algorithm, we can * just reject this value. This function then models a Zipf * distribution truncated to sys.maxint. */ - if (X > LONG_MAX) { - X = 0.0; /* X < 1 will be rejected */ + if (X > LONG_MAX || X < 1.0) { continue; } - if (X >= 1) { - T = pow(1.0 + 1.0/X, am1); + + T = pow(1.0 + 1.0/X, am1); + if (V*X*(T - 1.0)/(b - 1.0) <= T/b) { + return (long)X; } - } while ((X < 1) || ((V*X*(T-1.0)/(b-1.0)) > (T/b))); - return (long)X; + } } long rk_geometric_search(rk_state *state, double p) |