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



Attachments

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):
https://bitbucket.org/eigen/eigen/commits/d60e23b8af1
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) {
  n++;
  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;

with:

int xrand(n) {
  if(n==1)
   return std::rand();
  else
   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:
https://bitbucket.org/eigen/eigen/commits/00d109dc8349/
https://bitbucket.org/eigen/eigen/commits/8eb116bb13f9/ (3.2)

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