Line data Source code
1 : /* boost random/mersenne_twister.hpp header file
2 : *
3 : * Copyright Jens Maurer 2000-2001
4 : * Copyright Steven Watanabe 2010
5 : * Distributed under the Boost Software License, Version 1.0. (See
6 : * accompanying file LICENSE_1_0.txt or copy at
7 : * http://www.boost.org/LICENSE_1_0.txt)
8 : *
9 : * See http://www.boost.org for most recent version including documentation.
10 : *
11 : * $Id$
12 : *
13 : * Revision history
14 : * 2013-10-14 fixed some warnings with Wshadow (mgaunard)
15 : * 2001-02-18 moved to individual header files
16 : */
17 :
18 : #ifndef BOOST_RANDOM_MERSENNE_TWISTER_HPP
19 : #define BOOST_RANDOM_MERSENNE_TWISTER_HPP
20 :
21 : #include <iosfwd>
22 : #include <istream>
23 : #include <stdexcept>
24 : #include <boost/config.hpp>
25 : #include <boost/cstdint.hpp>
26 : #include <boost/integer/integer_mask.hpp>
27 : #include <boost/random/detail/config.hpp>
28 : #include <boost/random/detail/ptr_helper.hpp>
29 : #include <boost/random/detail/seed.hpp>
30 : #include <boost/random/detail/seed_impl.hpp>
31 : #include <boost/random/detail/generator_seed_seq.hpp>
32 : #include <boost/random/detail/polynomial.hpp>
33 :
34 : #include <boost/random/detail/disable_warnings.hpp>
35 :
36 : namespace boost {
37 : namespace random {
38 :
39 : /**
40 : * Instantiations of class template mersenne_twister_engine model a
41 : * \pseudo_random_number_generator. It uses the algorithm described in
42 : *
43 : * @blockquote
44 : * "Mersenne Twister: A 623-dimensionally equidistributed uniform
45 : * pseudo-random number generator", Makoto Matsumoto and Takuji Nishimura,
46 : * ACM Transactions on Modeling and Computer Simulation: Special Issue on
47 : * Uniform Random Number Generation, Vol. 8, No. 1, January 1998, pp. 3-30.
48 : * @endblockquote
49 : *
50 : * @xmlnote
51 : * The boost variant has been implemented from scratch and does not
52 : * derive from or use mt19937.c provided on the above WWW site. However, it
53 : * was verified that both produce identical output.
54 : * @endxmlnote
55 : *
56 : * The seeding from an integer was changed in April 2005 to address a
57 : * <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html">weakness</a>.
58 : *
59 : * The quality of the generator crucially depends on the choice of the
60 : * parameters. User code should employ one of the sensibly parameterized
61 : * generators such as \mt19937 instead.
62 : *
63 : * The generator requires considerable amounts of memory for the storage of
64 : * its state array. For example, \mt11213b requires about 1408 bytes and
65 : * \mt19937 requires about 2496 bytes.
66 : */
67 : template<class UIntType,
68 : std::size_t w, std::size_t n, std::size_t m, std::size_t r,
69 : UIntType a, std::size_t u, UIntType d, std::size_t s,
70 : UIntType b, std::size_t t,
71 : UIntType c, std::size_t l, UIntType f>
72 : class mersenne_twister_engine
73 : {
74 : public:
75 : typedef UIntType result_type;
76 : BOOST_STATIC_CONSTANT(std::size_t, word_size = w);
77 : BOOST_STATIC_CONSTANT(std::size_t, state_size = n);
78 : BOOST_STATIC_CONSTANT(std::size_t, shift_size = m);
79 : BOOST_STATIC_CONSTANT(std::size_t, mask_bits = r);
80 : BOOST_STATIC_CONSTANT(UIntType, xor_mask = a);
81 : BOOST_STATIC_CONSTANT(std::size_t, tempering_u = u);
82 : BOOST_STATIC_CONSTANT(UIntType, tempering_d = d);
83 : BOOST_STATIC_CONSTANT(std::size_t, tempering_s = s);
84 : BOOST_STATIC_CONSTANT(UIntType, tempering_b = b);
85 : BOOST_STATIC_CONSTANT(std::size_t, tempering_t = t);
86 : BOOST_STATIC_CONSTANT(UIntType, tempering_c = c);
87 : BOOST_STATIC_CONSTANT(std::size_t, tempering_l = l);
88 : BOOST_STATIC_CONSTANT(UIntType, initialization_multiplier = f);
89 : BOOST_STATIC_CONSTANT(UIntType, default_seed = 5489u);
90 :
91 : // backwards compatibility
92 : BOOST_STATIC_CONSTANT(UIntType, parameter_a = a);
93 : BOOST_STATIC_CONSTANT(std::size_t, output_u = u);
94 : BOOST_STATIC_CONSTANT(std::size_t, output_s = s);
95 : BOOST_STATIC_CONSTANT(UIntType, output_b = b);
96 : BOOST_STATIC_CONSTANT(std::size_t, output_t = t);
97 : BOOST_STATIC_CONSTANT(UIntType, output_c = c);
98 : BOOST_STATIC_CONSTANT(std::size_t, output_l = l);
99 :
100 : // old Boost.Random concept requirements
101 : BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
102 :
103 :
104 : /**
105 : * Constructs a @c mersenne_twister_engine and calls @c seed().
106 : */
107 0 : mersenne_twister_engine() { seed(); }
108 :
109 : /**
110 : * Constructs a @c mersenne_twister_engine and calls @c seed(value).
111 : */
112 : BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(mersenne_twister_engine,
113 : UIntType, value)
114 : { seed(value); }
115 : template<class It> mersenne_twister_engine(It& first, It last)
116 : { seed(first,last); }
117 :
118 : /**
119 : * Constructs a mersenne_twister_engine and calls @c seed(gen).
120 : *
121 : * @xmlnote
122 : * The copy constructor will always be preferred over
123 : * the templated constructor.
124 : * @endxmlnote
125 : */
126 : BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(mersenne_twister_engine,
127 : SeedSeq, seq)
128 : { seed(seq); }
129 :
130 : // compiler-generated copy ctor and assignment operator are fine
131 :
132 : /** Calls @c seed(default_seed). */
133 0 : void seed() { seed(default_seed); }
134 :
135 : /**
136 : * Sets the state x(0) to v mod 2w. Then, iteratively,
137 : * sets x(i) to
138 : * (i + f * (x(i-1) xor (x(i-1) rshift w-2))) mod 2<sup>w</sup>
139 : * for i = 1 .. n-1. x(n) is the first value to be returned by operator().
140 : */
141 : BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(mersenne_twister_engine, UIntType, value)
142 : {
143 : // New seeding algorithm from
144 : // http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html
145 : // In the previous versions, MSBs of the seed affected only MSBs of the
146 : // state x[].
147 : const UIntType mask = (max)();
148 : x[0] = value & mask;
149 : for (i = 1; i < n; i++) {
150 : // See Knuth "The Art of Computer Programming"
151 : // Vol. 2, 3rd ed., page 106
152 : x[i] = (f * (x[i-1] ^ (x[i-1] >> (w-2))) + i) & mask;
153 : }
154 :
155 : normalize_state();
156 : }
157 :
158 : /**
159 : * Seeds a mersenne_twister_engine using values produced by seq.generate().
160 : */
161 : BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(mersenne_twister_engine, SeeqSeq, seq)
162 : {
163 : detail::seed_array_int<w>(seq, x);
164 : i = n;
165 :
166 : normalize_state();
167 : }
168 :
169 : /** Sets the state of the generator using values from an iterator range. */
170 : template<class It>
171 : void seed(It& first, It last)
172 : {
173 : detail::fill_array_int<w>(first, last, x);
174 : i = n;
175 :
176 : normalize_state();
177 : }
178 :
179 : /** Returns the smallest value that the generator can produce. */
180 0 : static result_type min BOOST_PREVENT_MACRO_SUBSTITUTION ()
181 : { return 0; }
182 : /** Returns the largest value that the generator can produce. */
183 0 : static result_type max BOOST_PREVENT_MACRO_SUBSTITUTION ()
184 : { return boost::low_bits_mask_t<w>::sig_bits; }
185 :
186 : /** Produces the next value of the generator. */
187 : result_type operator()();
188 :
189 : /** Fills a range with random values */
190 : template<class Iter>
191 : void generate(Iter first, Iter last)
192 : { detail::generate_from_int(*this, first, last); }
193 :
194 : /**
195 : * Advances the state of the generator by @c z steps. Equivalent to
196 : *
197 : * @code
198 : * for(unsigned long long i = 0; i < z; ++i) {
199 : * gen();
200 : * }
201 : * @endcode
202 : */
203 : void discard(boost::uintmax_t z)
204 : {
205 : #ifndef BOOST_RANDOM_MERSENNE_TWISTER_DISCARD_THRESHOLD
206 : #define BOOST_RANDOM_MERSENNE_TWISTER_DISCARD_THRESHOLD 10000000
207 : #endif
208 : if(z > BOOST_RANDOM_MERSENNE_TWISTER_DISCARD_THRESHOLD) {
209 : discard_many(z);
210 : } else {
211 : for(boost::uintmax_t j = 0; j < z; ++j) {
212 : (*this)();
213 : }
214 : }
215 : }
216 :
217 : #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
218 : /** Writes a mersenne_twister_engine to a @c std::ostream */
219 : template<class CharT, class Traits>
220 : friend std::basic_ostream<CharT,Traits>&
221 : operator<<(std::basic_ostream<CharT,Traits>& os,
222 : const mersenne_twister_engine& mt)
223 : {
224 : mt.print(os);
225 : return os;
226 : }
227 :
228 : /** Reads a mersenne_twister_engine from a @c std::istream */
229 : template<class CharT, class Traits>
230 : friend std::basic_istream<CharT,Traits>&
231 : operator>>(std::basic_istream<CharT,Traits>& is,
232 : mersenne_twister_engine& mt)
233 : {
234 : for(std::size_t j = 0; j < mt.state_size; ++j)
235 : is >> mt.x[j] >> std::ws;
236 : // MSVC (up to 7.1) and Borland (up to 5.64) don't handle the template
237 : // value parameter "n" available from the class template scope, so use
238 : // the static constant with the same value
239 : mt.i = mt.state_size;
240 : return is;
241 : }
242 : #endif
243 :
244 : /**
245 : * Returns true if the two generators are in the same state,
246 : * and will thus produce identical sequences.
247 : */
248 : friend bool operator==(const mersenne_twister_engine& x_,
249 : const mersenne_twister_engine& y_)
250 : {
251 : if(x_.i < y_.i) return x_.equal_imp(y_);
252 : else return y_.equal_imp(x_);
253 : }
254 :
255 : /**
256 : * Returns true if the two generators are in different states.
257 : */
258 : friend bool operator!=(const mersenne_twister_engine& x_,
259 : const mersenne_twister_engine& y_)
260 : { return !(x_ == y_); }
261 :
262 : private:
263 : /// \cond show_private
264 :
265 : void twist();
266 :
267 : /**
268 : * Does the work of operator==. This is in a member function
269 : * for portability. Some compilers, such as msvc 7.1 and
270 : * Sun CC 5.10 can't access template parameters or static
271 : * members of the class from inline friend functions.
272 : *
273 : * requires i <= other.i
274 : */
275 : bool equal_imp(const mersenne_twister_engine& other) const
276 : {
277 : UIntType back[n];
278 : std::size_t offset = other.i - i;
279 : for(std::size_t j = 0; j + offset < n; ++j)
280 : if(x[j] != other.x[j+offset])
281 : return false;
282 : rewind(&back[n-1], offset);
283 : for(std::size_t j = 0; j < offset; ++j)
284 : if(back[j + n - offset] != other.x[j])
285 : return false;
286 : return true;
287 : }
288 :
289 : /**
290 : * Does the work of operator<<. This is in a member function
291 : * for portability.
292 : */
293 : template<class CharT, class Traits>
294 : void print(std::basic_ostream<CharT, Traits>& os) const
295 : {
296 : UIntType data[n];
297 : for(std::size_t j = 0; j < i; ++j) {
298 : data[j + n - i] = x[j];
299 : }
300 : if(i != n) {
301 : rewind(&data[n - i - 1], n - i);
302 : }
303 : os << data[0];
304 : for(std::size_t j = 1; j < n; ++j) {
305 : os << ' ' << data[j];
306 : }
307 : }
308 :
309 : /**
310 : * Copies z elements of the state preceding x[0] into
311 : * the array whose last element is last.
312 : */
313 : void rewind(UIntType* last, std::size_t z) const
314 : {
315 : const UIntType upper_mask = (~static_cast<UIntType>(0)) << r;
316 : const UIntType lower_mask = ~upper_mask;
317 : UIntType y0 = x[m-1] ^ x[n-1];
318 : if(y0 & (static_cast<UIntType>(1) << (w-1))) {
319 : y0 = ((y0 ^ a) << 1) | 1;
320 : } else {
321 : y0 = y0 << 1;
322 : }
323 : for(std::size_t sz = 0; sz < z; ++sz) {
324 : UIntType y1 =
325 : rewind_find(last, sz, m-1) ^ rewind_find(last, sz, n-1);
326 : if(y1 & (static_cast<UIntType>(1) << (w-1))) {
327 : y1 = ((y1 ^ a) << 1) | 1;
328 : } else {
329 : y1 = y1 << 1;
330 : }
331 : *(last - sz) = (y0 & upper_mask) | (y1 & lower_mask);
332 : y0 = y1;
333 : }
334 : }
335 :
336 : /**
337 : * Converts an arbitrary array into a valid generator state.
338 : * First we normalize x[0], so that it contains the same
339 : * value we would get by running the generator forwards
340 : * and then in reverse. (The low order r bits are redundant).
341 : * Then, if the state consists of all zeros, we set the
342 : * high order bit of x[0] to 1. This function only needs to
343 : * be called by seed, since the state transform preserves
344 : * this relationship.
345 : */
346 0 : void normalize_state()
347 : {
348 0 : const UIntType upper_mask = (~static_cast<UIntType>(0)) << r;
349 0 : const UIntType lower_mask = ~upper_mask;
350 0 : UIntType y0 = x[m-1] ^ x[n-1];
351 0 : if(y0 & (static_cast<UIntType>(1) << (w-1))) {
352 0 : y0 = ((y0 ^ a) << 1) | 1;
353 : } else {
354 0 : y0 = y0 << 1;
355 : }
356 0 : x[0] = (x[0] & upper_mask) | (y0 & lower_mask);
357 :
358 : // fix up the state if it's all zeroes.
359 0 : for(std::size_t j = 0; j < n; ++j) {
360 0 : if(x[j] != 0) return;
361 : }
362 0 : x[0] = static_cast<UIntType>(1) << (w-1);
363 : }
364 :
365 : /**
366 : * Given a pointer to the last element of the rewind array,
367 : * and the current size of the rewind array, finds an element
368 : * relative to the next available slot in the rewind array.
369 : */
370 : UIntType
371 : rewind_find(UIntType* last, std::size_t size, std::size_t j) const
372 : {
373 : std::size_t index = (j + n - size + n - 1) % n;
374 : if(index < n - size) {
375 : return x[index];
376 : } else {
377 : return *(last - (n - 1 - index));
378 : }
379 : }
380 :
381 : /**
382 : * Optimized algorithm for large jumps.
383 : *
384 : * Hiroshi Haramoto, Makoto Matsumoto, and Pierre L'Ecuyer. 2008.
385 : * A Fast Jump Ahead Algorithm for Linear Recurrences in a Polynomial
386 : * Space. In Proceedings of the 5th international conference on
387 : * Sequences and Their Applications (SETA '08).
388 : * DOI=10.1007/978-3-540-85912-3_26
389 : */
390 : void discard_many(boost::uintmax_t z)
391 : {
392 : // Compute the minimal polynomial, phi(t)
393 : // This depends only on the transition function,
394 : // which is constant. The characteristic
395 : // polynomial is the same as the minimal
396 : // polynomial for a maximum period generator
397 : // (which should be all specializations of
398 : // mersenne_twister.) Even if it weren't,
399 : // the characteristic polynomial is guaranteed
400 : // to be a multiple of the minimal polynomial,
401 : // which is good enough.
402 : detail::polynomial phi = get_characteristic_polynomial();
403 :
404 : // calculate g(t) = t^z % phi(t)
405 : detail::polynomial g = mod_pow_x(z, phi);
406 :
407 : // h(s_0, t) = \sum_{i=0}^{2k-1}o(s_i)t^{2k-i-1}
408 : detail::polynomial h;
409 : const std::size_t num_bits = w*n - r;
410 : for(std::size_t j = 0; j < num_bits * 2; ++j) {
411 : // Yes, we're advancing the generator state
412 : // here, but it doesn't matter because
413 : // we're going to overwrite it completely
414 : // in reconstruct_state.
415 : if(i >= n) twist();
416 : h[2*num_bits - j - 1] = x[i++] & UIntType(1);
417 : }
418 : // g(t)h(s_0, t)
419 : detail::polynomial gh = g * h;
420 : detail::polynomial result;
421 : for(std::size_t j = 0; j <= num_bits; ++j) {
422 : result[j] = gh[2*num_bits - j - 1];
423 : }
424 : reconstruct_state(result);
425 : }
426 : static detail::polynomial get_characteristic_polynomial()
427 : {
428 : const std::size_t num_bits = w*n - r;
429 : detail::polynomial helper;
430 : helper[num_bits - 1] = 1;
431 : mersenne_twister_engine tmp;
432 : tmp.reconstruct_state(helper);
433 : // Skip the first num_bits elements, since we
434 : // already know what they are.
435 : for(std::size_t j = 0; j < num_bits; ++j) {
436 : if(tmp.i >= n) tmp.twist();
437 : if(j == num_bits - 1)
438 : assert((tmp.x[tmp.i] & 1) == 1);
439 : else
440 : assert((tmp.x[tmp.i] & 1) == 0);
441 : ++tmp.i;
442 : }
443 : detail::polynomial phi;
444 : phi[num_bits] = 1;
445 : detail::polynomial next_bits = tmp.as_polynomial(num_bits);
446 : for(std::size_t j = 0; j < num_bits; ++j) {
447 : int val = next_bits[j] ^ phi[num_bits-j-1];
448 : phi[num_bits-j-1] = val;
449 : if(val) {
450 : for(std::size_t k = j + 1; k < num_bits; ++k) {
451 : phi[num_bits-k-1] ^= next_bits[k-j-1];
452 : }
453 : }
454 : }
455 : return phi;
456 : }
457 : detail::polynomial as_polynomial(std::size_t size) {
458 : detail::polynomial result;
459 : for(std::size_t j = 0; j < size; ++j) {
460 : if(i >= n) twist();
461 : result[j] = x[i++] & UIntType(1);
462 : }
463 : return result;
464 : }
465 : void reconstruct_state(const detail::polynomial& p)
466 : {
467 : const UIntType upper_mask = (~static_cast<UIntType>(0)) << r;
468 : const UIntType lower_mask = ~upper_mask;
469 : const std::size_t num_bits = w*n - r;
470 : for(std::size_t j = num_bits - n + 1; j <= num_bits; ++j)
471 : x[j % n] = p[j];
472 :
473 : UIntType y0 = 0;
474 : for(std::size_t j = num_bits + 1; j >= n - 1; --j) {
475 : UIntType y1 = x[j % n] ^ x[(j + m) % n];
476 : if(p[j - n + 1])
477 : y1 = (y1 ^ a) << UIntType(1) | UIntType(1);
478 : else
479 : y1 = y1 << UIntType(1);
480 : x[(j + 1) % n] = (y0 & upper_mask) | (y1 & lower_mask);
481 : y0 = y1;
482 : }
483 : i = 0;
484 : }
485 :
486 : /// \endcond
487 :
488 : // state representation: next output is o(x(i))
489 : // x[0] ... x[k] x[k+1] ... x[n-1] represents
490 : // x(i-k) ... x(i) x(i+1) ... x(i-k+n-1)
491 :
492 : UIntType x[n];
493 : std::size_t i;
494 : };
495 :
496 : /// \cond show_private
497 :
498 : #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
499 : // A definition is required even for integral static constants
500 : #define BOOST_RANDOM_MT_DEFINE_CONSTANT(type, name) \
501 : template<class UIntType, std::size_t w, std::size_t n, std::size_t m, \
502 : std::size_t r, UIntType a, std::size_t u, UIntType d, std::size_t s, \
503 : UIntType b, std::size_t t, UIntType c, std::size_t l, UIntType f> \
504 : const type mersenne_twister_engine<UIntType,w,n,m,r,a,u,d,s,b,t,c,l,f>::name
505 : BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, word_size);
506 : BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, state_size);
507 : BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, shift_size);
508 : BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, mask_bits);
509 : BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, xor_mask);
510 : BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, tempering_u);
511 : BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, tempering_d);
512 : BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, tempering_s);
513 : BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, tempering_b);
514 : BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, tempering_t);
515 : BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, tempering_c);
516 : BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, tempering_l);
517 : BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, initialization_multiplier);
518 : BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, default_seed);
519 : BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, parameter_a);
520 : BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, output_u );
521 : BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, output_s);
522 : BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, output_b);
523 : BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, output_t);
524 : BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, output_c);
525 : BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, output_l);
526 : BOOST_RANDOM_MT_DEFINE_CONSTANT(bool, has_fixed_range);
527 : #undef BOOST_RANDOM_MT_DEFINE_CONSTANT
528 : #endif
529 :
530 : template<class UIntType,
531 : std::size_t w, std::size_t n, std::size_t m, std::size_t r,
532 : UIntType a, std::size_t u, UIntType d, std::size_t s,
533 : UIntType b, std::size_t t,
534 : UIntType c, std::size_t l, UIntType f>
535 : void
536 0 : mersenne_twister_engine<UIntType,w,n,m,r,a,u,d,s,b,t,c,l,f>::twist()
537 : {
538 0 : const UIntType upper_mask = (~static_cast<UIntType>(0)) << r;
539 0 : const UIntType lower_mask = ~upper_mask;
540 :
541 0 : const std::size_t unroll_factor = 6;
542 0 : const std::size_t unroll_extra1 = (n-m) % unroll_factor;
543 0 : const std::size_t unroll_extra2 = (m-1) % unroll_factor;
544 :
545 : // split loop to avoid costly modulo operations
546 : { // extra scope for MSVC brokenness w.r.t. for scope
547 0 : for(std::size_t j = 0; j < n-m-unroll_extra1; j++) {
548 0 : UIntType y = (x[j] & upper_mask) | (x[j+1] & lower_mask);
549 0 : x[j] = x[j+m] ^ (y >> 1) ^ ((x[j+1]&1) * a);
550 : }
551 : }
552 : {
553 0 : for(std::size_t j = n-m-unroll_extra1; j < n-m; j++) {
554 0 : UIntType y = (x[j] & upper_mask) | (x[j+1] & lower_mask);
555 0 : x[j] = x[j+m] ^ (y >> 1) ^ ((x[j+1]&1) * a);
556 : }
557 : }
558 : {
559 0 : for(std::size_t j = n-m; j < n-1-unroll_extra2; j++) {
560 0 : UIntType y = (x[j] & upper_mask) | (x[j+1] & lower_mask);
561 0 : x[j] = x[j-(n-m)] ^ (y >> 1) ^ ((x[j+1]&1) * a);
562 : }
563 : }
564 : {
565 0 : for(std::size_t j = n-1-unroll_extra2; j < n-1; j++) {
566 : UIntType y = (x[j] & upper_mask) | (x[j+1] & lower_mask);
567 : x[j] = x[j-(n-m)] ^ (y >> 1) ^ ((x[j+1]&1) * a);
568 : }
569 : }
570 : // last iteration
571 0 : UIntType y = (x[n-1] & upper_mask) | (x[0] & lower_mask);
572 0 : x[n-1] = x[m-1] ^ (y >> 1) ^ ((x[0]&1) * a);
573 0 : i = 0;
574 0 : }
575 : /// \endcond
576 :
577 : template<class UIntType,
578 : std::size_t w, std::size_t n, std::size_t m, std::size_t r,
579 : UIntType a, std::size_t u, UIntType d, std::size_t s,
580 : UIntType b, std::size_t t,
581 : UIntType c, std::size_t l, UIntType f>
582 : inline typename
583 : mersenne_twister_engine<UIntType,w,n,m,r,a,u,d,s,b,t,c,l,f>::result_type
584 0 : mersenne_twister_engine<UIntType,w,n,m,r,a,u,d,s,b,t,c,l,f>::operator()()
585 : {
586 0 : if(i == n)
587 0 : twist();
588 : // Step 4
589 0 : UIntType z = x[i];
590 0 : ++i;
591 0 : z ^= ((z >> u) & d);
592 0 : z ^= ((z << s) & b);
593 0 : z ^= ((z << t) & c);
594 0 : z ^= (z >> l);
595 0 : return z;
596 : }
597 :
598 : /**
599 : * The specializations \mt11213b and \mt19937 are from
600 : *
601 : * @blockquote
602 : * "Mersenne Twister: A 623-dimensionally equidistributed
603 : * uniform pseudo-random number generator", Makoto Matsumoto
604 : * and Takuji Nishimura, ACM Transactions on Modeling and
605 : * Computer Simulation: Special Issue on Uniform Random Number
606 : * Generation, Vol. 8, No. 1, January 1998, pp. 3-30.
607 : * @endblockquote
608 : */
609 : typedef mersenne_twister_engine<uint32_t,32,351,175,19,0xccab8ee7,
610 : 11,0xffffffff,7,0x31b6ab00,15,0xffe50000,17,1812433253> mt11213b;
611 :
612 : /**
613 : * The specializations \mt11213b and \mt19937 are from
614 : *
615 : * @blockquote
616 : * "Mersenne Twister: A 623-dimensionally equidistributed
617 : * uniform pseudo-random number generator", Makoto Matsumoto
618 : * and Takuji Nishimura, ACM Transactions on Modeling and
619 : * Computer Simulation: Special Issue on Uniform Random Number
620 : * Generation, Vol. 8, No. 1, January 1998, pp. 3-30.
621 : * @endblockquote
622 : */
623 : typedef mersenne_twister_engine<uint32_t,32,624,397,31,0x9908b0df,
624 : 11,0xffffffff,7,0x9d2c5680,15,0xefc60000,18,1812433253> mt19937;
625 :
626 : #if !defined(BOOST_NO_INT64_T) && !defined(BOOST_NO_INTEGRAL_INT64_T)
627 : typedef mersenne_twister_engine<uint64_t,64,312,156,31,
628 : UINT64_C(0xb5026f5aa96619e9),29,UINT64_C(0x5555555555555555),17,
629 : UINT64_C(0x71d67fffeda60000),37,UINT64_C(0xfff7eee000000000),43,
630 : UINT64_C(6364136223846793005)> mt19937_64;
631 : #endif
632 :
633 : /// \cond show_deprecated
634 :
635 : template<class UIntType,
636 : int w, int n, int m, int r,
637 : UIntType a, int u, std::size_t s,
638 : UIntType b, int t,
639 : UIntType c, int l, UIntType v>
640 : class mersenne_twister :
641 : public mersenne_twister_engine<UIntType,
642 : w, n, m, r, a, u, ~(UIntType)0, s, b, t, c, l, 1812433253>
643 : {
644 : typedef mersenne_twister_engine<UIntType,
645 : w, n, m, r, a, u, ~(UIntType)0, s, b, t, c, l, 1812433253> base_type;
646 : public:
647 : mersenne_twister() {}
648 : BOOST_RANDOM_DETAIL_GENERATOR_CONSTRUCTOR(mersenne_twister, Gen, gen)
649 : { seed(gen); }
650 : BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(mersenne_twister, UIntType, val)
651 : { seed(val); }
652 : template<class It>
653 : mersenne_twister(It& first, It last) : base_type(first, last) {}
654 : void seed() { base_type::seed(); }
655 : BOOST_RANDOM_DETAIL_GENERATOR_SEED(mersenne_twister, Gen, gen)
656 : {
657 : detail::generator_seed_seq<Gen> seq(gen);
658 : base_type::seed(seq);
659 : }
660 : BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(mersenne_twister, UIntType, val)
661 : { base_type::seed(val); }
662 : template<class It>
663 : void seed(It& first, It last) { base_type::seed(first, last); }
664 : };
665 :
666 : /// \endcond
667 :
668 : } // namespace random
669 :
670 : using random::mt11213b;
671 : using random::mt19937;
672 : using random::mt19937_64;
673 :
674 : } // namespace boost
675 :
676 : BOOST_RANDOM_PTR_HELPER_SPEC(boost::mt11213b)
677 : BOOST_RANDOM_PTR_HELPER_SPEC(boost::mt19937)
678 : BOOST_RANDOM_PTR_HELPER_SPEC(boost::mt19937_64)
679 :
680 : #include <boost/random/detail/enable_warnings.hpp>
681 :
682 : #endif // BOOST_RANDOM_MERSENNE_TWISTER_HPP
|
