GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-07-04 04:05:19 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fix possible overflow and biais in integer random generator

Browse Source
This commit is contained in:
Gael Guennebaud 2016-08-23 13:25:31 +02:00
parent 581b6472d1
commit 82147cefff
2 changed files with 42 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
Eigen/src/Core/MathFunctions.h
View File

@ -615,16 +615,18 @@ struct random_default_impl<Scalar, false, true>
typedef typename conditional<NumTraits<Scalar>::IsSigned,std::ptrdiff_t,std::size_t>::type ScalarX; typedef typename conditional<NumTraits<Scalar>::IsSigned,std::ptrdiff_t,std::size_t>::type ScalarX;
if(y<x) if(y<x)
return x; return x;
// the following difference might overflow on a 32 bits system,
// but since y>=x the result converted to an unsigned long is still correct.
std::size_t range = ScalarX(y)-ScalarX(x); std::size_t range = ScalarX(y)-ScalarX(x);
std::size_t offset = 0; std::size_t offset = 0;
// rejection sampling // rejection sampling
std::size_t divisor = (range+RAND_MAX-1)/(range+1); std::size_t divisor = 1;
std::size_t multiplier = (range+RAND_MAX-1)/std::size_t(RAND_MAX); std::size_t multiplier = 1;
if(range<RAND_MAX) divisor = (std::size_t(RAND_MAX)+1)/(range+1);
else multiplier = 1 + range/(std::size_t(RAND_MAX)+1);
do { do {
offset = ( (std::size_t(std::rand()) * multiplier) / divisor ); offset = (std::size_t(std::rand()) * multiplier) / divisor;
} while (offset > range); } while (offset > range);
return Scalar(ScalarX(x) + offset); return Scalar(ScalarX(x) + offset);
} }

29
test/rand.cpp
View File

@ -9,6 +9,8 @@
#include "main.h" #include "main.h"
typedef long long int64;
template<typename Scalar> Scalar check_in_range(Scalar x, Scalar y) template<typename Scalar> Scalar check_in_range(Scalar x, Scalar y)
{ {
Scalar r = internal::random<Scalar>(x,y); Scalar r = internal::random<Scalar>(x,y);
@ -35,13 +37,30 @@ template<typename Scalar> void check_all_in_range(Scalar x, Scalar y)
VERIFY( (mask>0).all() ); VERIFY( (mask>0).all() );
} }
template<typename Scalar> void check_histogram(Scalar x, Scalar y, int bins)
{
Array<int,1,Dynamic> hist(bins);
hist.fill(0);
int f = 100000;
int n = bins*f;
int64 range = int64(y)-int64(x);
int divisor = int((range+1)/bins);
assert(((range+1)%bins)==0);
for(int k=0; k<n; ++k)
{
Scalar r = check_in_range(x,y);
hist( int((int64(r)-int64(x))/divisor) )++;
}
VERIFY( (((hist.cast<double>()/double(f))-1.0).abs()<0.02).all() );
}
void test_rand() void test_rand()
{ {
long long_ref = NumTraits<long>::highest()/10; long long_ref = NumTraits<long>::highest()/10;
signed char char_offset = (std::min)(g_repeat,64); signed char char_offset = (std::min)(g_repeat,64);
signed char short_offset = (std::min)(g_repeat,16000); signed char short_offset = (std::min)(g_repeat,16000);
for(int i = 0; i < g_repeat*10; i++) { for(int i = 0; i < g_repeat*10000; i++) {
CALL_SUBTEST(check_in_range<float>(10,11)); CALL_SUBTEST(check_in_range<float>(10,11));
CALL_SUBTEST(check_in_range<float>(1.24234523,1.24234523)); CALL_SUBTEST(check_in_range<float>(1.24234523,1.24234523));
CALL_SUBTEST(check_in_range<float>(-1,1)); CALL_SUBTEST(check_in_range<float>(-1,1));
@ -56,6 +75,7 @@ void test_rand()
CALL_SUBTEST(check_in_range<short>(0,-1)); CALL_SUBTEST(check_in_range<short>(0,-1));
CALL_SUBTEST(check_in_range<long>(0,-1)); CALL_SUBTEST(check_in_range<long>(0,-1));
CALL_SUBTEST(check_in_range<int>(-673456,673456)); CALL_SUBTEST(check_in_range<int>(-673456,673456));
CALL_SUBTEST(check_in_range<int>(-RAND_MAX+10,RAND_MAX-10));
CALL_SUBTEST(check_in_range<short>(-24345,24345)); CALL_SUBTEST(check_in_range<short>(-24345,24345));
CALL_SUBTEST(check_in_range<long>(-long_ref,long_ref)); CALL_SUBTEST(check_in_range<long>(-long_ref,long_ref));
} }
@ -88,4 +108,11 @@ void test_rand()
CALL_SUBTEST(check_all_in_range<long>(-11-g_repeat,-11)); CALL_SUBTEST(check_all_in_range<long>(-11-g_repeat,-11));
CALL_SUBTEST(check_all_in_range<long>(-long_ref,-long_ref+g_repeat)); CALL_SUBTEST(check_all_in_range<long>(-long_ref,-long_ref+g_repeat));
CALL_SUBTEST(check_all_in_range<long>( long_ref, long_ref+g_repeat)); CALL_SUBTEST(check_all_in_range<long>( long_ref, long_ref+g_repeat));
CALL_SUBTEST(check_histogram<int>(-5,5,11));
int bins = 100;
CALL_SUBTEST(check_histogram<int>(-3333,-3333+bins*(3333/bins)-1,bins));
bins = 1000;
CALL_SUBTEST(check_histogram<int>(-RAND_MAX+10,-RAND_MAX+10+bins*(RAND_MAX/bins)-1,bins));
CALL_SUBTEST(check_histogram<int>(-RAND_MAX+10,-int64(RAND_MAX)+10+bins*(2*int64(RAND_MAX)/bins)-1,bins));
} }