New user self-registration is currently disabled. Please email eigen-core-team @ if you need an account.
Bug 1275 - upper bound breach in random<short>
upper bound breach in random<short>
Product: Eigen
Classification: Unclassified
Component: Core - general
All All
: Normal Unknown
Assigned To: Nobody
Depends on:
  Show dependency treegraph
Reported: 2016-08-15 12:04 UTC by Ilya Mezhirov
Modified: 2016-08-23 11:27 UTC (History)
3 users (show)


Description Ilya Mezhirov 2016-08-15 12:04:15 UTC
Sometimes random<short>(0, n) returns n + 1. This could be responsible for some out of bounds errors in the tests (I've tracked a failure in sparse_basic_1 to this).

The problem is in this formula (Eigen/src/Core/MathFunctions.h):

  static inline Scalar run(const Scalar& x, const Scalar& y)
    return x + Scalar((NonInteger(y)-x+1) * std::rand() / (RAND_MAX + NonInteger(1)));

When Scalar is a 16-bit int, NonInteger is a 32-bit float with 24 bits for the mantissa. However, std::rand() is wider than Scalar since it produces 31 meaningful bits. So when std::rand() returns a value that is very close to RAND_MAX, 24-bit mantissa cannot distinguish them.

Here's a plain C snippet that shows this formula in action: 

#include <stdio.h>
int main() {
    float f = (float) (18 + 1) * 2147483640 / (2147483647 + (float) 1);
    printf("%f %d\n", f, (int) f);
    return 0;

which prints "19.000000 19".

The most obvious solution is to replace NonInteger with double for this case. But there's probably some similar effect for 64-bit integers and doubles. So perhaps you might want to treat 64-bit Scalar as a special case and replace NonInteger with uint64_t.
Comment 1 Christoph Hertzberg 2016-08-15 13:06:31 UTC
I copied the improved rand_impl from the devel branch (as well as the corresponding unit test):
Comment 2 Ilya Mezhirov 2016-08-15 13:59:26 UTC
That was quick, thanks!
Please note that 64 bit size of size_t is not guaranteed though.
Comment 3 Gael Guennebaud 2016-08-22 12:44:29 UTC
Indeed, another approach would be to recursively call rand() when range>RAND_MAX:

int n = 1;
int rand_max = RAND_MAX;
while(rand_max<range) {
  rand_max = rand_max*(RAND_MAX+1)+RAND_MAX;
divisor = ((rand_max+1)/range);
do { x = xrand(n); } while(x>=range * divisor);
return x / divisor;


int xrand(n) {
   return std::rand();
   return xrand(n-1) * (RAND_MAX+1) + std::rand();
Comment 4 Gael Guennebaud 2016-08-23 11:27:50 UTC
Regarding rand, I fixed some overflow issues: (3.2)

Note You need to log in before you can comment on or make changes to this bug.