#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;
-
- uint64_t cache; /* assumed to be twice as long as the usual random number */
- int bits_left; /* bit results */
-
+ std::tr1::mt19937 mt;
+ uint64_t supremum;
+ uint64_t RAND_SUP;
+ double supremum_reciprocal;
+ double RAND_SUP_REC;
+
+ // 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
#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
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.
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. */
operator()(uint64_t local_supremum)
{
return scale_down(mt(), local_supremum,
- (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;
+ double(local_supremum * RAND_SUP_REC));
}
- /** @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
bits_left -= bits;
if (bits_left < 32)
{
- cache |= (((uint64_t)mt()) << bits_left);
+ cache |= ((uint64_t(mt())) << bits_left);
bits_left += 32;
}
return res;
OSDN Git Service
