New user self-registration is disabled due to spam. Please email eigen-core-team @ lists.tuxfamily.org if you need an account.
Before reporting a bug, please make sure that your Eigen version is up-to-date!
Bug 1275 - upper bound breach in random<short>
upper bound breach in random<short>
 Status: RESOLVED FIXED None Eigen Unclassified Core - general (show other bugs) 3.2 All All Normal Unknown Nobody

 Reported: 2016-08-15 12:04 UTC by Ilya Mezhirov 2016-08-23 11:27 UTC (History) 3 users (show) chtz gael.guennebaud jacob.benoit.1

Attachments

 Ilya Mezhirov 2016-08-15 12:04:15 UTC ```Sometimes random(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 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.``` 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``` 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.``` 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 * divisor); return x / divisor; with: int xrand(n) { if(n==1) return std::rand(); else return xrand(n-1) * (RAND_MAX+1) + std::rand(); }``` 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.