[Numpy-svn] r5364 - branches/1.1.x/numpy/random/mtrand

[numpy-svn at scipy.org]

Author: rkern
Date: 2008-07-08 20:56:22 -0500 (Tue, 08 Jul 2008)
New Revision: 5364

Modified:
   branches/1.1.x/numpy/random/mtrand/distributions.c
Log:
BUG: Make sure the Zipf results are within the allowable range.

Modified: branches/1.1.x/numpy/random/mtrand/distributions.c
===================================================================
--- branches/1.1.x/numpy/random/mtrand/distributions.c	2008-07-09 01:55:33 UTC (rev 5363)
+++ branches/1.1.x/numpy/random/mtrand/distributions.c	2008-07-09 01:56:22 UTC (rev 5364)
@@ -676,16 +676,23 @@
 {
     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;
 }