This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen
Bug 1275 - upper bound breach in random<short>
Summary: upper bound breach in random<short>
Status: RESOLVED FIXED
Alias: None
Product: Eigen
Classification: Unclassified
Component: Core - general (show other bugs)
Version: 3.2
Hardware: All All
: Normal Unknown
Assignee: Nobody
URL:
Whiteboard:
Keywords:
Depends on:
Blocks:
 
Reported: 2016-08-15 12:04 UTC by Ilya Mezhirov
Modified: 2019-12-04 16:06 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)
Comment 5 Nobody 2019-12-04 16:06:17 UTC
-- GitLab Migration Automatic Message --

This bug has been migrated to gitlab.com's GitLab instance and has been closed from further activity.

You can subscribe and participate further through the new bug through this link to our GitLab instance: https://gitlab.com/libeigen/eigen/issues/1275.

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