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.

1 files 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()