3 // Copyright (C) 2007, 2008 Free Software Foundation, Inc.
5 // This file is part of the GNU ISO C++ Library. This library is free
6 // software; you can redistribute it and/or modify it under the terms
7 // of the GNU General Public License as published by the Free Software
8 // Foundation; either version 2, or (at your option) any later
11 // This library is distributed in the hope that it will be useful, but
12 // WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // General Public License for more details.
16 // You should have received a copy of the GNU General Public License
17 // along with this library; see the file COPYING. If not, write to
18 // the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
19 // MA 02111-1307, USA.
21 // As a special exception, you may use this file as part of a free
22 // software library without restriction. Specifically, if other files
23 // instantiate templates or use macros or inline functions from this
24 // file, or you compile this file and link it with other files to
25 // produce an executable, this file does not by itself cause the
26 // resulting executable to be covered by the GNU General Public
27 // License. This exception does not however invalidate any other
28 // reasons why the executable file might be covered by the GNU General
31 /** @file parallel/random_number.h
32 * @brief Random number generator based on the Mersenne twister.
33 * This file is a GNU parallel extension to the Standard C++ Library.
36 // Written by Johannes Singler.
38 #ifndef _GLIBCXX_PARALLEL_RANDOM_NUMBER_H
39 #define _GLIBCXX_PARALLEL_RANDOM_NUMBER_H 1
41 #include <parallel/types.h>
43 namespace __gnu_parallel
45 // XXX use tr1 random number.
46 // http://www.math.keio.ac.jp/matumoto/emt.html
47 template<typename UIntType
, int w
, int n
, int m
, int r
, UIntType a
, int u
,
48 int s
, UIntType b
, int t
, UIntType c
, int l
, UIntType val
>
49 class mersenne_twister
52 typedef UIntType result_type
;
53 static const int word_size
= w
;
54 static const int state_size
= n
;
55 static const int shift_size
= m
;
56 static const int mask_bits
= r
;
57 static const UIntType parameter_a
= a
;
58 static const int output_u
= u
;
59 static const int output_s
= s
;
60 static const UIntType output_b
= b
;
61 static const int output_t
= t
;
62 static const UIntType output_c
= c
;
63 static const int output_l
= l
;
65 static const bool has_fixed_range
= false;
67 mersenne_twister() { seed(); }
69 #if defined(__SUNPRO_CC) && (__SUNPRO_CC <= 0x520)
70 // Work around overload resolution problem (Gennadiy E. Rozental)
72 mersenne_twister(const UIntType
& value
)
75 mersenne_twister(UIntType value
)
80 mersenne_twister(It
& first
, It last
)
83 template<typename Generator
>
85 mersenne_twister(Generator
& gen
)
88 // compiler-generated copy ctor and assignment operator are fine
92 { seed(UIntType(5489)); }
94 #if defined(__SUNPRO_CC) && (__SUNPRO_CC <= 0x520)
95 // Work around overload resolution problem (Gennadiy E. Rozental)
97 seed(const UIntType
& value
)
103 // New seeding algorithm from
104 // http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html
105 // In the previous versions, MSBs of the seed affected only MSBs of the
107 const UIntType mask
= ~0u;
109 for (i
= 1; i
< n
; ++i
)
111 // See Knuth "The Art of Computer Programming" Vol. 2,
113 x
[i
] = (1812433253UL * (x
[i
-1] ^ (x
[i
-1] >> (w
-2))) + i
) & mask
;
117 // For GCC, moving this function out-of-line prevents inlining, which may
118 // reduce overall object code size. However, MSVC does not grok
119 // out-of-line definitions of member function templates.
120 template<typename Generator
>
122 seed(Generator
& gen
)
124 // I could have used std::generate_n, but it takes "gen" by value
125 for (int j
= 0; j
< n
; ++j
)
130 template<typename It
>
132 seed(It
& first
, It last
)
135 for (j
= 0; j
< n
&& first
!= last
; ++j
, ++first
)
138 /* if (first == last && j < n)
139 throw std::invalid_argument("mersenne_twister::seed");*/
149 // avoid "left shift count >= with of type" warning
151 for (int i
= 0; i
< w
; ++i
)
160 validation(result_type v
)
163 #ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
166 operator==(const mersenne_twister
& x
, const mersenne_twister
& y
)
168 for (int j
= 0; j
< state_size
; ++j
)
169 if (x
.compute(j
) != y
.compute(j
))
175 operator!=(const mersenne_twister
& x
, const mersenne_twister
& y
)
176 { return !(x
== y
); }
178 // Use a member function; Streamable concept not supported.
180 operator==(const mersenne_twister
& rhs
) const
182 for (int j
= 0; j
< state_size
; ++j
)
183 if (compute(j
) != rhs
.compute(j
))
189 operator!=(const mersenne_twister
& rhs
) const
190 { return !(*this == rhs
); }
194 // returns x(i-n+index), where index is in 0..n-1
196 compute(unsigned int index
) const
198 // equivalent to (i-n+index) % 2n, but doesn't produce negative numbers
199 return x
[ (i
+ n
+ index
) % (2*n
) ];
205 // state representation: next output is o(x(i))
206 // x[0] ... x[k] x[k+1] ... x[n-1] x[n] ... x[2*n-1] represents
207 // x(i-k) ... x(i) x(i+1) ... x(i-k+n-1) x(i-k-n) ... x[i(i-k-1)]
208 // The goal is to always have x(i-n) ... x(i-1) available for
209 // operator== and save/restore.
215 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
216 // A definition is required even for integral static constants
217 template<typename UIntType
, int w
, int n
, int m
, int r
, UIntType a
, int u
,
218 int s
, UIntType b
, int t
, UIntType c
, int l
, UIntType val
>
220 mersenne_twister
<UIntType
,w
,n
,m
,r
,a
,u
,s
,b
,t
,c
,l
,val
>::has_fixed_range
;
222 template<typename UIntType
, int w
, int n
, int m
, int r
, UIntType a
, int u
,
223 int s
, UIntType b
, int t
, UIntType c
, int l
, UIntType val
>
225 mersenne_twister
<UIntType
,w
,n
,m
,r
,a
,u
,s
,b
,t
,c
,l
,val
>::state_size
;
227 template<typename UIntType
, int w
, int n
, int m
, int r
, UIntType a
, int u
,
228 int s
, UIntType b
, int t
, UIntType c
, int l
, UIntType val
>
230 mersenne_twister
<UIntType
,w
,n
,m
,r
,a
,u
,s
,b
,t
,c
,l
,val
>::shift_size
;
232 template<typename UIntType
, int w
, int n
, int m
, int r
, UIntType a
, int u
,
233 int s
, UIntType b
, int t
, UIntType c
, int l
, UIntType val
>
235 mersenne_twister
<UIntType
,w
,n
,m
,r
,a
,u
,s
,b
,t
,c
,l
,val
>::mask_bits
;
237 template<typename UIntType
, int w
, int n
, int m
, int r
, UIntType a
, int u
,
238 int s
, UIntType b
, int t
, UIntType c
, int l
, UIntType val
>
240 mersenne_twister
<UIntType
,w
,n
,m
,r
,a
,u
,s
,b
,t
,c
,l
,val
>::parameter_a
;
242 template<typename UIntType
, int w
, int n
, int m
, int r
, UIntType a
, int u
,
243 int s
, UIntType b
, int t
, UIntType c
, int l
, UIntType val
>
245 mersenne_twister
<UIntType
,w
,n
,m
,r
,a
,u
,s
,b
,t
,c
,l
,val
>::output_u
;
247 template<typename UIntType
, int w
, int n
, int m
, int r
, UIntType a
, int u
,
248 int s
, UIntType b
, int t
, UIntType c
, int l
, UIntType val
>
250 mersenne_twister
<UIntType
,w
,n
,m
,r
,a
,u
,s
,b
,t
,c
,l
,val
>::output_s
;
252 template<typename UIntType
, int w
, int n
, int m
, int r
, UIntType a
, int u
,
253 int s
, UIntType b
, int t
, UIntType c
, int l
, UIntType val
>
255 mersenne_twister
<UIntType
,w
,n
,m
,r
,a
,u
,s
,b
,t
,c
,l
,val
>::output_b
;
257 template<typename UIntType
, int w
, int n
, int m
, int r
, UIntType a
, int u
,
258 int s
, UIntType b
, int t
, UIntType c
, int l
, UIntType val
>
260 mersenne_twister
<UIntType
,w
,n
,m
,r
,a
,u
,s
,b
,t
,c
,l
,val
>::output_t
;
262 template<typename UIntType
, int w
, int n
, int m
, int r
, UIntType a
, int u
,
263 int s
, UIntType b
, int t
, UIntType c
, int l
, UIntType val
>
265 mersenne_twister
<UIntType
,w
,n
,m
,r
,a
,u
,s
,b
,t
,c
,l
,val
>::output_c
;
267 template<typename UIntType
, int w
, int n
, int m
, int r
, UIntType a
, int u
,
268 int s
, UIntType b
, int t
, UIntType c
, int l
, UIntType val
>
270 mersenne_twister
<UIntType
,w
,n
,m
,r
,a
,u
,s
,b
,t
,c
,l
,val
>::output_l
;
273 template<typename UIntType
, int w
, int n
, int m
, int r
, UIntType a
, int u
,
274 int s
, UIntType b
, int t
, UIntType c
, int l
, UIntType val
>
276 mersenne_twister
<UIntType
,w
,n
,m
,r
,a
,u
,s
,b
,t
,c
,l
,val
>::twist(int block
)
278 const UIntType upper_mask
= (~0u) << r
;
279 const UIntType lower_mask
= ~upper_mask
;
283 for (int j
= n
; j
< 2*n
; ++j
)
285 UIntType y
= (x
[j
-n
] & upper_mask
) | (x
[j
-(n
-1)] & lower_mask
);
286 x
[j
] = x
[j
-(n
-m
)] ^ (y
>> 1) ^ (y
&1 ? a
: 0);
291 // split loop to avoid costly modulo operations
292 { // extra scope for MSVC brokenness w.r.t. for scope
293 for (int j
= 0; j
< n
-m
; ++j
)
295 UIntType y
= (x
[j
+n
] & upper_mask
) | (x
[j
+n
+1] & lower_mask
);
296 x
[j
] = x
[j
+n
+m
] ^ (y
>> 1) ^ (y
&1 ? a
: 0);
300 for (int j
= n
-m
; j
< n
-1; ++j
)
302 UIntType y
= (x
[j
+n
] & upper_mask
) | (x
[j
+n
+1] & lower_mask
);
303 x
[j
] = x
[j
-(n
-m
)] ^ (y
>> 1) ^ (y
&1 ? a
: 0);
306 UIntType y
= (x
[2*n
-1] & upper_mask
) | (x
[0] & lower_mask
);
307 x
[n
-1] = x
[m
-1] ^ (y
>> 1) ^ (y
&1 ? a
: 0);
312 template<typename UIntType
, int w
, int n
, int m
, int r
, UIntType a
, int u
,
313 int s
, UIntType b
, int t
, UIntType c
, int l
, UIntType val
>
315 typename mersenne_twister
<UIntType
,w
,n
,m
,r
,a
,u
,s
,b
,t
,c
,l
,val
>::result_type
316 mersenne_twister
<UIntType
,w
,n
,m
,r
,a
,u
,s
,b
,t
,c
,l
,val
>::operator()()
333 typedef mersenne_twister
<uint32_t,32,351,175,19,0xccab8ee7,11,
334 7,0x31b6ab00,15,0xffe50000,17, 0xa37d3c92> mt11213b
;
336 // validation by experiment from mt19937.c
337 typedef mersenne_twister
<uint32_t,32,624,397,31,0x9908b0df,11,
338 7,0x9d2c5680,15,0xefc60000,18, 3346425566U> mt19937
;
340 /** @brief Random number generator, based on the Mersenne twister. */
345 uint64_t supremum
, RAND_SUP
;
346 double supremum_reciprocal
, RAND_SUP_REC
;
348 uint64_t cache
; /* assumed to be twice as long as the usual random number */
349 int bits_left
; /* bit results */
352 scale_down(uint64_t x
,
353 #if _GLIBCXX_SCALE_DOWN_FPU
354 uint64_t /*supremum*/, double supremum_reciprocal
)
356 uint64_t supremum
, double /*supremum_reciprocal*/)
359 #if _GLIBCXX_SCALE_DOWN_FPU
360 return (uint32_t)(x
* supremum_reciprocal
);
362 return static_cast<uint32_t>(x
% supremum
);
367 /** @brief Default constructor. Seed with 0. */
369 : mt(0), supremum(0x100000000ULL
),
370 RAND_SUP(1ULL << (sizeof(uint32_t) * 8)),
371 supremum_reciprocal((double)supremum
/ (double)RAND_SUP
),
372 RAND_SUP_REC(1.0 / (double)RAND_SUP
),
373 cache(0), bits_left(0) { }
375 /** @brief Constructor.
376 * @param seed Random seed.
377 * @param supremum Generate integer random numbers in the
378 * interval @c [0,supremum). */
379 random_number(uint32_t seed
, uint64_t supremum
= 0x100000000ULL
)
380 : mt(seed
), supremum(supremum
),
381 RAND_SUP(1ULL << (sizeof(uint32_t) * 8)),
382 supremum_reciprocal((double)supremum
/ (double)RAND_SUP
),
383 RAND_SUP_REC(1.0 / (double)RAND_SUP
),
384 cache(0), bits_left(0) { }
386 /** @brief Generate unsigned random 32-bit integer. */
389 { return scale_down(mt(), supremum
, supremum_reciprocal
); }
391 /** @brief Generate unsigned random 32-bit integer in the
392 interval @c [0,local_supremum). */
394 operator()(uint64_t local_supremum
)
396 return scale_down(mt(), local_supremum
,
397 (double)local_supremum
* RAND_SUP_REC
);
400 /** @brief Set the random seed.
401 * @param seed to set. */
403 set_seed(uint32_t seed
)
410 /** @brief Generate a number of random bits, compile-time parameter. */
415 unsigned long res
= cache
& ((1 << bits
) - 1);
416 cache
= cache
>> bits
;
420 cache
|= (((uint64_t)mt()) << bits_left
);
426 /** @brief Generate a number of random bits, run-time parameter.
427 * @param bits Number of bits to generate. */
429 genrand_bits(int bits
)
431 unsigned long res
= cache
& ((1 << bits
) - 1);
432 cache
= cache
>> bits
;
436 cache
|= (((uint64_t)mt()) << bits_left
);
443 } // namespace __gnu_parallel