From d99ab35f9e886a014be6d47606d232af1e668f76 Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Fri, 13 Mar 2015 21:12:46 +0100
Subject: [PATCH] Fix internal::random(x,y) for integer types. The previous
 implementation could return y+1. The new implementation uses rejection
 sampling to get an unbiased behabior.

---
 Eigen/src/Core/MathFunctions.h | 21 +++++++++++++++++++--
 1 file changed, 19 insertions(+), 2 deletions(-)

diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h
index 878f38e92..3c76a58b9 100644
--- a/Eigen/src/Core/MathFunctions.h
+++ b/Eigen/src/Core/MathFunctions.h
@@ -525,8 +525,25 @@ struct random_default_impl<Scalar, false, true>
   typedef typename NumTraits<Scalar>::NonInteger NonInteger;
 
   static inline Scalar run(const Scalar& x, const Scalar& y)
-  {
-    return x + Scalar((NonInteger(y)-x+1) * std::rand() / (RAND_MAX + NonInteger(1)));
+  { 
+    using std::max;
+    Scalar range = (max)(Scalar(0),Scalar(y-x));
+    Scalar offset = 0;
+    if(range<=RAND_MAX)
+    {
+      // rejection sampling
+      int divisor = RAND_MAX/(range+1);
+
+      do {
+        offset = Scalar(std::rand() / divisor);
+      } while (offset > range);
+    }
+    else
+    {
+      offset = std::rand() * range;
+    }
+    
+    return x + offset;
   }
 
   static inline Scalar run()