From 8989a31d6ddaaa1f7df0e5da9bdee88cd4900a66 Mon Sep 17 00:00:00 2001
From: Yang Li <yangli2@gmail.com>
Date: Tue, 12 Sep 2017 10:01:48 -0400
Subject: BUG: Fix possibly undefined cast of double -> long.

Current double to long casting in the zipf function depends on
non-standardized behavior when the double is too big to fit in a long.
This is potentially dangerous and makes the code fail with tools such as
AddressSanitizer.

Checks are added here to prevent overflow during casting and make sure
we get the desired behavior.
---
 numpy/random/mtrand/distributions.c | 24 +++++++++++++++---------
 1 file changed, 15 insertions(+), 9 deletions(-)

diff --git a/numpy/random/mtrand/distributions.c b/numpy/random/mtrand/distributions.c
index 7673f92b4..a455140f8 100644
--- a/numpy/random/mtrand/distributions.c
+++ b/numpy/random/mtrand/distributions.c
@@ -45,6 +45,7 @@
 #include <stdio.h>
 #include <math.h>
 #include <stdlib.h>
+#include <limits.h>
 
 #ifndef min
 #define min(x,y) ((x<y)?x:y)
@@ -719,26 +720,31 @@ double rk_wald(rk_state *state, double mean, double scale)
 
 long rk_zipf(rk_state *state, double a)
 {
-    double T, U, V;
-    long X;
+    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);
         V = rk_double(state);
-        X = (long)floor(pow(U, -1.0/am1));
+        X = 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
+         * 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.
          */
-        T = pow(1.0 + 1.0/X, am1);
-    } while (((V*X*(T-1.0)/(b-1.0)) > (T/b)) || X < 1);
-    return X;
+        if (X > LONG_MAX) {
+            X = 0.0;    /* X < 1 will be rejected */
+            continue;
+        }
+        if (X >= 1) {
+            T = pow(1.0 + 1.0/X, am1);
+        }
+    } 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)
-- 
cgit v1.2.3