Fix internal::random(x,y) for integer types. The previous implementation could return y+1. The new implementation uses rejection sampling to get an unbiased behabior.

This commit is contained in:
Gael Guennebaud 2015-03-13 21:12:46 +01:00
parent 8580eb6808
commit d99ab35f9e

View File

@ -525,8 +525,25 @@ struct random_default_impl<Scalar, false, true>
typedef typename NumTraits<Scalar>::NonInteger NonInteger;
static inline Scalar run(const Scalar& x, const Scalar& y)
{
return x + Scalar((NonInteger(y)-x+1) * std::rand() / (RAND_MAX + NonInteger(1)));
{
using std::max;
Scalar range = (max)(Scalar(0),Scalar(y-x));
Scalar offset = 0;
if(range<=RAND_MAX)
{
// rejection sampling
int divisor = RAND_MAX/(range+1);
do {
offset = Scalar(std::rand() / divisor);
} while (offset > range);
}
else
{
offset = std::rand() * range;
}
return x + offset;
}
static inline Scalar run()