random_number.h: Use TR1's mersenne_twister.
2008-02-16 Benjamin Kosnik <bkoz@redhat.com> * include/parallel/random_number.h: Use TR1's mersenne_twister. (random_number::genrand_bits()): Remove. (random_number::set_seed): Remove. From-SVN: r132371
This commit is contained in:
Benjamin Kosnik
Benjamin Kosnik
parent
9e2d69391e
commit
459af5a02d
@ -1,3 +1,9 @@
|
||||
2008-02-16 Benjamin Kosnik <bkoz@redhat.com>
|
||||
|
||||
* include/parallel/random_number.h: Use TR1's mersenne_twister.
|
||||
(random_number::genrand_bits()): Remove.
|
||||
(random_number::set_seed): Remove.
|
||||
|
||||
2008-02-15 Benjamin Kosnik <bkoz@redhat.com>
|
||||
|
||||
* include/parallel/types.h: Remove enum parallelism.
|
||||
|
||||
@ -39,315 +39,26 @@
|
||||
#define _GLIBCXX_PARALLEL_RANDOM_NUMBER_H 1
|
||||
|
||||
#include <parallel/types.h>
|
||||
#include <tr1/random>
|
||||
|
||||
namespace __gnu_parallel
|
||||
{
|
||||
// XXX use tr1 random number.
|
||||
// http://www.math.keio.ac.jp/matumoto/emt.html
|
||||
template<typename UIntType, int w, int n, int m, int r, UIntType a, int u,
|
||||
int s, UIntType b, int t, UIntType c, int l, UIntType val>
|
||||
class mersenne_twister
|
||||
{
|
||||
public:
|
||||
typedef UIntType result_type;
|
||||
static const int word_size = w;
|
||||
static const int state_size = n;
|
||||
static const int shift_size = m;
|
||||
static const int mask_bits = r;
|
||||
static const UIntType parameter_a = a;
|
||||
static const int output_u = u;
|
||||
static const int output_s = s;
|
||||
static const UIntType output_b = b;
|
||||
static const int output_t = t;
|
||||
static const UIntType output_c = c;
|
||||
static const int output_l = l;
|
||||
|
||||
static const bool has_fixed_range = false;
|
||||
|
||||
mersenne_twister() { seed(); }
|
||||
|
||||
#if defined(__SUNPRO_CC) && (__SUNPRO_CC <= 0x520)
|
||||
// Work around overload resolution problem (Gennadiy E. Rozental)
|
||||
explicit
|
||||
mersenne_twister(const UIntType& value)
|
||||
#else
|
||||
explicit
|
||||
mersenne_twister(UIntType value)
|
||||
#endif
|
||||
{ seed(value); }
|
||||
|
||||
template<typename It>
|
||||
mersenne_twister(It& first, It last)
|
||||
{ seed(first,last); }
|
||||
|
||||
template<typename Generator>
|
||||
explicit
|
||||
mersenne_twister(Generator & gen)
|
||||
{ seed(gen); }
|
||||
|
||||
// compiler-generated copy ctor and assignment operator are fine
|
||||
|
||||
void
|
||||
seed()
|
||||
{ seed(UIntType(5489)); }
|
||||
|
||||
#if defined(__SUNPRO_CC) && (__SUNPRO_CC <= 0x520)
|
||||
// Work around overload resolution problem (Gennadiy E. Rozental)
|
||||
void
|
||||
seed(const UIntType& value)
|
||||
#else
|
||||
void
|
||||
seed(UIntType value)
|
||||
#endif
|
||||
{
|
||||
// New seeding algorithm from
|
||||
// http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html
|
||||
// In the previous versions, MSBs of the seed affected only MSBs of the
|
||||
// state x[].
|
||||
const UIntType mask = ~0u;
|
||||
x[0] = value & mask;
|
||||
for (i = 1; i < n; ++i)
|
||||
{
|
||||
// See Knuth "The Art of Computer Programming" Vol. 2,
|
||||
// 3rd ed., page 106
|
||||
x[i] = (1812433253UL * (x[i-1] ^ (x[i-1] >> (w-2))) + i) & mask;
|
||||
}
|
||||
}
|
||||
|
||||
// For GCC, moving this function out-of-line prevents inlining, which may
|
||||
// reduce overall object code size. However, MSVC does not grok
|
||||
// out-of-line definitions of member function templates.
|
||||
template<typename Generator>
|
||||
void
|
||||
seed(Generator & gen)
|
||||
{
|
||||
// I could have used std::generate_n, but it takes "gen" by value
|
||||
for (int j = 0; j < n; ++j)
|
||||
x[j] = gen();
|
||||
i = n;
|
||||
}
|
||||
|
||||
template<typename It>
|
||||
void
|
||||
seed(It& first, It last)
|
||||
{
|
||||
int j;
|
||||
for (j = 0; j < n && first != last; ++j, ++first)
|
||||
x[j] = *first;
|
||||
i = n;
|
||||
/* if (first == last && j < n)
|
||||
throw std::invalid_argument("mersenne_twister::seed");*/
|
||||
}
|
||||
|
||||
result_type
|
||||
min() const
|
||||
{ return 0; }
|
||||
|
||||
result_type
|
||||
max() const
|
||||
{
|
||||
// avoid "left shift count >= with of type" warning
|
||||
result_type res = 0;
|
||||
for (int i = 0; i < w; ++i)
|
||||
res |= (1u << i);
|
||||
return res;
|
||||
}
|
||||
|
||||
result_type
|
||||
operator()();
|
||||
|
||||
static bool
|
||||
validation(result_type v)
|
||||
{ return val == v; }
|
||||
|
||||
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
|
||||
|
||||
friend bool
|
||||
operator==(const mersenne_twister& x, const mersenne_twister& y)
|
||||
{
|
||||
for (int j = 0; j < state_size; ++j)
|
||||
if (x.compute(j) != y.compute(j))
|
||||
return false;
|
||||
return true;
|
||||
}
|
||||
|
||||
friend bool
|
||||
operator!=(const mersenne_twister& x, const mersenne_twister& y)
|
||||
{ return !(x == y); }
|
||||
#else
|
||||
// Use a member function; Streamable concept not supported.
|
||||
bool
|
||||
operator==(const mersenne_twister& rhs) const
|
||||
{
|
||||
for (int j = 0; j < state_size; ++j)
|
||||
if (compute(j) != rhs.compute(j))
|
||||
return false;
|
||||
return true;
|
||||
}
|
||||
|
||||
bool
|
||||
operator!=(const mersenne_twister& rhs) const
|
||||
{ return !(*this == rhs); }
|
||||
#endif
|
||||
|
||||
private:
|
||||
// returns x(i-n+index), where index is in 0..n-1
|
||||
UIntType
|
||||
compute(unsigned int index) const
|
||||
{
|
||||
// equivalent to (i-n+index) % 2n, but doesn't produce negative numbers
|
||||
return x[ (i + n + index) % (2*n) ];
|
||||
}
|
||||
|
||||
void
|
||||
twist(int block);
|
||||
|
||||
// state representation: next output is o(x(i))
|
||||
// x[0] ... x[k] x[k+1] ... x[n-1] x[n] ... x[2*n-1] represents
|
||||
// x(i-k) ... x(i) x(i+1) ... x(i-k+n-1) x(i-k-n) ... x[i(i-k-1)]
|
||||
// The goal is to always have x(i-n) ... x(i-1) available for
|
||||
// operator== and save/restore.
|
||||
|
||||
UIntType x[2*n];
|
||||
int i;
|
||||
};
|
||||
|
||||
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
|
||||
// A definition is required even for integral static constants
|
||||
template<typename UIntType, int w, int n, int m, int r, UIntType a, int u,
|
||||
int s, UIntType b, int t, UIntType c, int l, UIntType val>
|
||||
const bool
|
||||
mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::has_fixed_range;
|
||||
|
||||
template<typename UIntType, int w, int n, int m, int r, UIntType a, int u,
|
||||
int s, UIntType b, int t, UIntType c, int l, UIntType val>
|
||||
const int
|
||||
mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::state_size;
|
||||
|
||||
template<typename UIntType, int w, int n, int m, int r, UIntType a, int u,
|
||||
int s, UIntType b, int t, UIntType c, int l, UIntType val>
|
||||
const int
|
||||
mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::shift_size;
|
||||
|
||||
template<typename UIntType, int w, int n, int m, int r, UIntType a, int u,
|
||||
int s, UIntType b, int t, UIntType c, int l, UIntType val>
|
||||
const int
|
||||
mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::mask_bits;
|
||||
|
||||
template<typename UIntType, int w, int n, int m, int r, UIntType a, int u,
|
||||
int s, UIntType b, int t, UIntType c, int l, UIntType val>
|
||||
const UIntType
|
||||
mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::parameter_a;
|
||||
|
||||
template<typename UIntType, int w, int n, int m, int r, UIntType a, int u,
|
||||
int s, UIntType b, int t, UIntType c, int l, UIntType val>
|
||||
const int
|
||||
mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_u;
|
||||
|
||||
template<typename UIntType, int w, int n, int m, int r, UIntType a, int u,
|
||||
int s, UIntType b, int t, UIntType c, int l, UIntType val>
|
||||
const int
|
||||
mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_s;
|
||||
|
||||
template<typename UIntType, int w, int n, int m, int r, UIntType a, int u,
|
||||
int s, UIntType b, int t, UIntType c, int l, UIntType val>
|
||||
const UIntType
|
||||
mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_b;
|
||||
|
||||
template<typename UIntType, int w, int n, int m, int r, UIntType a, int u,
|
||||
int s, UIntType b, int t, UIntType c, int l, UIntType val>
|
||||
const int
|
||||
mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_t;
|
||||
|
||||
template<typename UIntType, int w, int n, int m, int r, UIntType a, int u,
|
||||
int s, UIntType b, int t, UIntType c, int l, UIntType val>
|
||||
const UIntType
|
||||
mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_c;
|
||||
|
||||
template<typename UIntType, int w, int n, int m, int r, UIntType a, int u,
|
||||
int s, UIntType b, int t, UIntType c, int l, UIntType val>
|
||||
const int
|
||||
mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_l;
|
||||
#endif
|
||||
|
||||
template<typename UIntType, int w, int n, int m, int r, UIntType a, int u,
|
||||
int s, UIntType b, int t, UIntType c, int l, UIntType val>
|
||||
void
|
||||
mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::twist(int block)
|
||||
{
|
||||
const UIntType upper_mask = (~0u) << r;
|
||||
const UIntType lower_mask = ~upper_mask;
|
||||
|
||||
if (block == 0)
|
||||
{
|
||||
for (int j = n; j < 2*n; ++j)
|
||||
{
|
||||
UIntType y = (x[j-n] & upper_mask) | (x[j-(n-1)] & lower_mask);
|
||||
x[j] = x[j-(n-m)] ^ (y >> 1) ^ (y&1 ? a : 0);
|
||||
}
|
||||
}
|
||||
else if (block == 1)
|
||||
{
|
||||
// split loop to avoid costly modulo operations
|
||||
{ // extra scope for MSVC brokenness w.r.t. for scope
|
||||
for (int j = 0; j < n-m; ++j)
|
||||
{
|
||||
UIntType y = (x[j+n] & upper_mask) | (x[j+n+1] & lower_mask);
|
||||
x[j] = x[j+n+m] ^ (y >> 1) ^ (y&1 ? a : 0);
|
||||
}
|
||||
}
|
||||
|
||||
for (int j = n-m; j < n-1; ++j)
|
||||
{
|
||||
UIntType y = (x[j+n] & upper_mask) | (x[j+n+1] & lower_mask);
|
||||
x[j] = x[j-(n-m)] ^ (y >> 1) ^ (y&1 ? a : 0);
|
||||
}
|
||||
// last iteration
|
||||
UIntType y = (x[2*n-1] & upper_mask) | (x[0] & lower_mask);
|
||||
x[n-1] = x[m-1] ^ (y >> 1) ^ (y&1 ? a : 0);
|
||||
i = 0;
|
||||
}
|
||||
}
|
||||
|
||||
template<typename UIntType, int w, int n, int m, int r, UIntType a, int u,
|
||||
int s, UIntType b, int t, UIntType c, int l, UIntType val>
|
||||
inline
|
||||
typename mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::result_type
|
||||
mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::operator()()
|
||||
{
|
||||
if (i == n)
|
||||
twist(0);
|
||||
else if (i >= 2*n)
|
||||
twist(1);
|
||||
// Step 4
|
||||
UIntType z = x[i];
|
||||
++i;
|
||||
z ^= (z >> u);
|
||||
z ^= ((z << s) & b);
|
||||
z ^= ((z << t) & c);
|
||||
z ^= (z >> l);
|
||||
return z;
|
||||
}
|
||||
|
||||
|
||||
typedef mersenne_twister<uint32_t,32,351,175,19,0xccab8ee7,11,
|
||||
7,0x31b6ab00,15,0xffe50000,17, 0xa37d3c92> mt11213b;
|
||||
|
||||
// validation by experiment from mt19937.c
|
||||
typedef mersenne_twister<uint32_t,32,624,397,31,0x9908b0df,11,
|
||||
7,0x9d2c5680,15,0xefc60000,18, 3346425566U> mt19937;
|
||||
|
||||
/** @brief Random number generator, based on the Mersenne twister. */
|
||||
class random_number
|
||||
{
|
||||
private:
|
||||
mt19937 mt;
|
||||
uint64_t supremum, RAND_SUP;
|
||||
double supremum_reciprocal, RAND_SUP_REC;
|
||||
std::tr1::mt19937 mt;
|
||||
uint64_t supremum;
|
||||
uint64_t RAND_SUP;
|
||||
double supremum_reciprocal;
|
||||
double RAND_SUP_REC;
|
||||
|
||||
uint64_t cache; /* assumed to be twice as long as the usual random number */
|
||||
int bits_left; /* bit results */
|
||||
// Assumed to be twice as long as the usual random number.
|
||||
uint64_t cache;
|
||||
|
||||
// Bit results.
|
||||
int bits_left;
|
||||
|
||||
static uint32_t
|
||||
scale_down(uint64_t x,
|
||||
#if _GLIBCXX_SCALE_DOWN_FPU
|
||||
@ -357,7 +68,7 @@ namespace __gnu_parallel
|
||||
#endif
|
||||
{
|
||||
#if _GLIBCXX_SCALE_DOWN_FPU
|
||||
return (uint32_t)(x * supremum_reciprocal);
|
||||
return uint32_t(x * supremum_reciprocal);
|
||||
#else
|
||||
return static_cast<uint32_t>(x % supremum);
|
||||
#endif
|
||||
@ -368,8 +79,8 @@ namespace __gnu_parallel
|
||||
random_number()
|
||||
: mt(0), supremum(0x100000000ULL),
|
||||
RAND_SUP(1ULL << (sizeof(uint32_t) * 8)),
|
||||
supremum_reciprocal((double)supremum / (double)RAND_SUP),
|
||||
RAND_SUP_REC(1.0 / (double)RAND_SUP),
|
||||
supremum_reciprocal(double(supremum) / double(RAND_SUP)),
|
||||
RAND_SUP_REC(1.0 / double(RAND_SUP)),
|
||||
cache(0), bits_left(0) { }
|
||||
|
||||
/** @brief Constructor.
|
||||
@ -379,8 +90,8 @@ namespace __gnu_parallel
|
||||
random_number(uint32_t seed, uint64_t supremum = 0x100000000ULL)
|
||||
: mt(seed), supremum(supremum),
|
||||
RAND_SUP(1ULL << (sizeof(uint32_t) * 8)),
|
||||
supremum_reciprocal((double)supremum / (double)RAND_SUP),
|
||||
RAND_SUP_REC(1.0 / (double)RAND_SUP),
|
||||
supremum_reciprocal(double(supremum) / double(RAND_SUP)),
|
||||
RAND_SUP_REC(1.0 / double(RAND_SUP)),
|
||||
cache(0), bits_left(0) { }
|
||||
|
||||
/** @brief Generate unsigned random 32-bit integer. */
|
||||
@ -394,35 +105,9 @@ namespace __gnu_parallel
|
||||
operator()(uint64_t local_supremum)
|
||||
{
|
||||
return scale_down(mt(), local_supremum,
|
||||
(double)local_supremum * RAND_SUP_REC);
|
||||
double(local_supremum * RAND_SUP_REC));
|
||||
}
|
||||
|
||||
/** @brief Set the random seed.
|
||||
* @param seed to set. */
|
||||
void
|
||||
set_seed(uint32_t seed)
|
||||
{
|
||||
mt.seed(seed);
|
||||
cache = mt();
|
||||
bits_left = 32;
|
||||
}
|
||||
|
||||
/** @brief Generate a number of random bits, compile-time parameter. */
|
||||
template<int bits>
|
||||
unsigned long
|
||||
genrand_bits()
|
||||
{
|
||||
unsigned long res = cache & ((1 << bits) - 1);
|
||||
cache = cache >> bits;
|
||||
bits_left -= bits;
|
||||
if (bits_left < 32)
|
||||
{
|
||||
cache |= (((uint64_t)mt()) << bits_left);
|
||||
bits_left += 32;
|
||||
}
|
||||
return res;
|
||||
}
|
||||
|
||||
/** @brief Generate a number of random bits, run-time parameter.
|
||||
* @param bits Number of bits to generate. */
|
||||
unsigned long
|
||||
@ -433,7 +118,7 @@ namespace __gnu_parallel
|
||||
bits_left -= bits;
|
||||
if (bits_left < 32)
|
||||
{
|
||||
cache |= (((uint64_t)mt()) << bits_left);
|
||||
cache |= ((uint64_t(mt())) << bits_left);
|
||||
bits_left += 32;
|
||||
}
|
||||
return res;
|
||||
|
||||
Loading…
Reference in New Issue
Block a user