* @brief Random number generator based on the Mersenne twister.
* This file is a GNU parallel extension to the Standard C++ Library.
*/
#ifndef _GLIBCXX_PARALLEL_RANDOM_NUMBER_H
#define _GLIBCXX_PARALLEL_RANDOM_NUMBER_H 1
#include <parallel/types.h>
#include <tr1/random>
namespace __gnu_parallel
{
class random_number
{
private:
std::tr1::mt19937 mt;
uint64 supremum;
uint64 RAND_SUP;
double supremum_reciprocal;
double RAND_SUP_REC;
uint64 cache;
int bits_left;
static uint32
scale_down(uint64 x,
#if _GLIBCXX_SCALE_DOWN_FPU
uint64 , double supremum_reciprocal)
#else
uint64 supremum, double )
#endif
{
#if _GLIBCXX_SCALE_DOWN_FPU
return uint32(x * supremum_reciprocal);
#else
return static_cast<uint32>(x % supremum);
#endif
}
public:
random_number()
: mt(0), supremum(0x100000000ULL),
RAND_SUP(1ULL << (sizeof(uint32) * 8)),
supremum_reciprocal(double(supremum) / double(RAND_SUP)),
RAND_SUP_REC(1.0 / double(RAND_SUP)),
cache(0), bits_left(0) { }
* @param seed Random seed.
* @param supremum Generate integer random numbers in the
* interval @c [0,supremum). */
random_number(uint32 seed, uint64 supremum = 0x100000000ULL)
: mt(seed), supremum(supremum),
RAND_SUP(1ULL << (sizeof(uint32) * 8)),
supremum_reciprocal(double(supremum) / double(RAND_SUP)),
RAND_SUP_REC(1.0 / double(RAND_SUP)),
cache(0), bits_left(0) { }
uint32
operator()()
{ return scale_down(mt(), supremum, supremum_reciprocal); }
interval @c [0,local_supremum). */
uint32
operator()(uint64 local_supremum)
{
return scale_down(mt(), local_supremum,
double(local_supremum * RAND_SUP_REC));
}
* @param bits Number of bits to generate. */
unsigned long
genrand_bits(int bits)
{
unsigned long res = cache & ((1 << bits) - 1);
cache = cache >> bits;
bits_left -= bits;
if (bits_left < 32)
{
cache |= ((uint64(mt())) << bits_left);
bits_left += 32;
}
return res;
}
};
}
#endif