BUG: Make sure the Zipf results are within the allowable range.

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;
 }