]> git.donarmstrong.com Git - rsem.git/blob - boost/random/mersenne_twister.hpp
Updated boost to v1.55.0
[rsem.git] / boost / random / mersenne_twister.hpp
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: mersenne_twister.hpp 74867 2011-10-09 23:13:31Z steven_watanabe $
12  *
13  * Revision history
14  *  2001-02-18  moved to individual header files
15  */
16
17 #ifndef BOOST_RANDOM_MERSENNE_TWISTER_HPP
18 #define BOOST_RANDOM_MERSENNE_TWISTER_HPP
19
20 #include <iosfwd>
21 #include <istream>
22 #include <stdexcept>
23 #include <boost/config.hpp>
24 #include <boost/cstdint.hpp>
25 #include <boost/integer/integer_mask.hpp>
26 #include <boost/random/detail/config.hpp>
27 #include <boost/random/detail/ptr_helper.hpp>
28 #include <boost/random/detail/seed.hpp>
29 #include <boost/random/detail/seed_impl.hpp>
30 #include <boost/random/detail/generator_seed_seq.hpp>
31
32 namespace boost {
33 namespace random {
34
35 /**
36  * Instantiations of class template mersenne_twister_engine model a
37  * \pseudo_random_number_generator. It uses the algorithm described in
38  *
39  *  @blockquote
40  *  "Mersenne Twister: A 623-dimensionally equidistributed uniform
41  *  pseudo-random number generator", Makoto Matsumoto and Takuji Nishimura,
42  *  ACM Transactions on Modeling and Computer Simulation: Special Issue on
43  *  Uniform Random Number Generation, Vol. 8, No. 1, January 1998, pp. 3-30. 
44  *  @endblockquote
45  *
46  * @xmlnote
47  * The boost variant has been implemented from scratch and does not
48  * derive from or use mt19937.c provided on the above WWW site. However, it
49  * was verified that both produce identical output.
50  * @endxmlnote
51  *
52  * The seeding from an integer was changed in April 2005 to address a
53  * <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html">weakness</a>.
54  * 
55  * The quality of the generator crucially depends on the choice of the
56  * parameters.  User code should employ one of the sensibly parameterized
57  * generators such as \mt19937 instead.
58  *
59  * The generator requires considerable amounts of memory for the storage of
60  * its state array. For example, \mt11213b requires about 1408 bytes and
61  * \mt19937 requires about 2496 bytes.
62  */
63 template<class UIntType,
64          std::size_t w, std::size_t n, std::size_t m, std::size_t r,
65          UIntType a, std::size_t u, UIntType d, std::size_t s,
66          UIntType b, std::size_t t,
67          UIntType c, std::size_t l, UIntType f>
68 class mersenne_twister_engine
69 {
70 public:
71     typedef UIntType result_type;
72     BOOST_STATIC_CONSTANT(std::size_t, word_size = w);
73     BOOST_STATIC_CONSTANT(std::size_t, state_size = n);
74     BOOST_STATIC_CONSTANT(std::size_t, shift_size = m);
75     BOOST_STATIC_CONSTANT(std::size_t, mask_bits = r);
76     BOOST_STATIC_CONSTANT(UIntType, xor_mask = a);
77     BOOST_STATIC_CONSTANT(std::size_t, tempering_u = u);
78     BOOST_STATIC_CONSTANT(UIntType, tempering_d = d);
79     BOOST_STATIC_CONSTANT(std::size_t, tempering_s = s);
80     BOOST_STATIC_CONSTANT(UIntType, tempering_b = b);
81     BOOST_STATIC_CONSTANT(std::size_t, tempering_t = t);
82     BOOST_STATIC_CONSTANT(UIntType, tempering_c = c);
83     BOOST_STATIC_CONSTANT(std::size_t, tempering_l = l);
84     BOOST_STATIC_CONSTANT(UIntType, initialization_multiplier = f);
85     BOOST_STATIC_CONSTANT(UIntType, default_seed = 5489u);
86   
87     // backwards compatibility
88     BOOST_STATIC_CONSTANT(UIntType, parameter_a = a);
89     BOOST_STATIC_CONSTANT(std::size_t, output_u = u);
90     BOOST_STATIC_CONSTANT(std::size_t, output_s = s);
91     BOOST_STATIC_CONSTANT(UIntType, output_b = b);
92     BOOST_STATIC_CONSTANT(std::size_t, output_t = t);
93     BOOST_STATIC_CONSTANT(UIntType, output_c = c);
94     BOOST_STATIC_CONSTANT(std::size_t, output_l = l);
95     
96     // old Boost.Random concept requirements
97     BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
98
99
100     /**
101      * Constructs a @c mersenne_twister_engine and calls @c seed().
102      */
103     mersenne_twister_engine() { seed(); }
104
105     /**
106      * Constructs a @c mersenne_twister_engine and calls @c seed(value).
107      */
108     BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(mersenne_twister_engine,
109                                                UIntType, value)
110     { seed(value); }
111     template<class It> mersenne_twister_engine(It& first, It last)
112     { seed(first,last); }
113
114     /**
115      * Constructs a mersenne_twister_engine and calls @c seed(gen).
116      *
117      * @xmlnote
118      * The copy constructor will always be preferred over
119      * the templated constructor.
120      * @endxmlnote
121      */
122     BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(mersenne_twister_engine,
123                                              SeedSeq, seq)
124     { seed(seq); }
125
126     // compiler-generated copy ctor and assignment operator are fine
127
128     /** Calls @c seed(default_seed). */
129     void seed() { seed(default_seed); }
130
131     /**
132      * Sets the state x(0) to v mod 2w. Then, iteratively,
133      * sets x(i) to
134      * (i + f * (x(i-1) xor (x(i-1) rshift w-2))) mod 2<sup>w</sup>
135      * for i = 1 .. n-1. x(n) is the first value to be returned by operator().
136      */
137     BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(mersenne_twister_engine, UIntType, value)
138     {
139         // New seeding algorithm from 
140         // http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html
141         // In the previous versions, MSBs of the seed affected only MSBs of the
142         // state x[].
143         const UIntType mask = (max)();
144         x[0] = value & mask;
145         for (i = 1; i < n; i++) {
146             // See Knuth "The Art of Computer Programming"
147             // Vol. 2, 3rd ed., page 106
148             x[i] = (f * (x[i-1] ^ (x[i-1] >> (w-2))) + i) & mask;
149         }
150     }
151     
152     /**
153      * Seeds a mersenne_twister_engine using values produced by seq.generate().
154      */
155     BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(mersenne_twister_engine, SeeqSeq, seq)
156     {
157         detail::seed_array_int<w>(seq, x);
158         i = n;
159
160         // fix up the state if it's all zeroes.
161         if((x[0] & (~static_cast<UIntType>(0) << r)) == 0) {
162             for(std::size_t j = 1; j < n; ++j) {
163                 if(x[j] != 0) return;
164             }
165             x[0] = static_cast<UIntType>(1) << (w-1);
166         }
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         // fix up the state if it's all zeroes.
177         if((x[0] & (~static_cast<UIntType>(0) << r)) == 0) {
178             for(std::size_t j = 1; j < n; ++j) {
179                 if(x[j] != 0) return;
180             }
181             x[0] = static_cast<UIntType>(1) << (w-1);
182         }
183     }
184   
185     /** Returns the smallest value that the generator can produce. */
186     static result_type min BOOST_PREVENT_MACRO_SUBSTITUTION ()
187     { return 0; }
188     /** Returns the largest value that the generator can produce. */
189     static result_type max BOOST_PREVENT_MACRO_SUBSTITUTION ()
190     { return boost::low_bits_mask_t<w>::sig_bits; }
191     
192     /** Produces the next value of the generator. */
193     result_type operator()();
194
195     /** Fills a range with random values */
196     template<class Iter>
197     void generate(Iter first, Iter last)
198     { detail::generate_from_int(*this, first, last); }
199
200     /**
201      * Advances the state of the generator by @c z steps.  Equivalent to
202      *
203      * @code
204      * for(unsigned long long i = 0; i < z; ++i) {
205      *     gen();
206      * }
207      * @endcode
208      */
209     void discard(boost::uintmax_t z)
210     {
211         for(boost::uintmax_t j = 0; j < z; ++j) {
212             (*this)();
213         }
214     }
215
216 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
217     /** Writes a mersenne_twister_engine to a @c std::ostream */
218     template<class CharT, class Traits>
219     friend std::basic_ostream<CharT,Traits>&
220     operator<<(std::basic_ostream<CharT,Traits>& os,
221                const mersenne_twister_engine& mt)
222     {
223         mt.print(os);
224         return os;
225     }
226     
227     /** Reads a mersenne_twister_engine from a @c std::istream */
228     template<class CharT, class Traits>
229     friend std::basic_istream<CharT,Traits>&
230     operator>>(std::basic_istream<CharT,Traits>& is,
231                mersenne_twister_engine& mt)
232     {
233         for(std::size_t j = 0; j < mt.state_size; ++j)
234             is >> mt.x[j] >> std::ws;
235         // MSVC (up to 7.1) and Borland (up to 5.64) don't handle the template
236         // value parameter "n" available from the class template scope, so use
237         // the static constant with the same value
238         mt.i = mt.state_size;
239         return is;
240     }
241 #endif
242
243     /**
244      * Returns true if the two generators are in the same state,
245      * and will thus produce identical sequences.
246      */
247     friend bool operator==(const mersenne_twister_engine& x,
248                            const mersenne_twister_engine& y)
249     {
250         if(x.i < y.i) return x.equal_imp(y);
251         else return y.equal_imp(x);
252     }
253     
254     /**
255      * Returns true if the two generators are in different states.
256      */
257     friend bool operator!=(const mersenne_twister_engine& x,
258                            const mersenne_twister_engine& y)
259     { return !(x == y); }
260
261 private:
262     /// \cond show_private
263
264     void twist();
265
266     /**
267      * Does the work of operator==.  This is in a member function
268      * for portability.  Some compilers, such as msvc 7.1 and
269      * Sun CC 5.10 can't access template parameters or static
270      * members of the class from inline friend functions.
271      *
272      * requires i <= other.i
273      */
274     bool equal_imp(const mersenne_twister_engine& other) const
275     {
276         UIntType back[n];
277         std::size_t offset = other.i - i;
278         for(std::size_t j = 0; j + offset < n; ++j)
279             if(x[j] != other.x[j+offset])
280                 return false;
281         rewind(&back[n-1], offset);
282         for(std::size_t j = 0; j < offset; ++j)
283             if(back[j + n - offset] != other.x[j])
284                 return false;
285         return true;
286     }
287
288     /**
289      * Does the work of operator<<.  This is in a member function
290      * for portability.
291      */
292     template<class CharT, class Traits>
293     void print(std::basic_ostream<CharT, Traits>& os) const
294     {
295         UIntType data[n];
296         for(std::size_t j = 0; j < i; ++j) {
297             data[j + n - i] = x[j];
298         }
299         if(i != n) {
300             rewind(&data[n - i - 1], n - i);
301         }
302         os << data[0];
303         for(std::size_t j = 1; j < n; ++j) {
304             os << ' ' << data[j];
305         }
306     }
307
308     /**
309      * Copies z elements of the state preceding x[0] into
310      * the array whose last element is last.
311      */
312     void rewind(UIntType* last, std::size_t z) const
313     {
314         const UIntType upper_mask = (~static_cast<UIntType>(0)) << r;
315         const UIntType lower_mask = ~upper_mask;
316         UIntType y0 = x[m-1] ^ x[n-1];
317         if(y0 & (static_cast<UIntType>(1) << (w-1))) {
318             y0 = ((y0 ^ a) << 1) | 1;
319         } else {
320             y0 = y0 << 1;
321         }
322         for(std::size_t sz = 0; sz < z; ++sz) {
323             UIntType y1 =
324                 rewind_find(last, sz, m-1) ^ rewind_find(last, sz, n-1);
325             if(y1 & (static_cast<UIntType>(1) << (w-1))) {
326                 y1 = ((y1 ^ a) << 1) | 1;
327             } else {
328                 y1 = y1 << 1;
329             }
330             *(last - sz) = (y0 & upper_mask) | (y1 & lower_mask);
331             y0 = y1;
332         }
333     }
334
335     /**
336      * Given a pointer to the last element of the rewind array,
337      * and the current size of the rewind array, finds an element
338      * relative to the next available slot in the rewind array.
339      */
340     UIntType
341     rewind_find(UIntType* last, std::size_t size, std::size_t j) const
342     {
343         std::size_t index = (j + n - size + n - 1) % n;
344         if(index < n - size) {
345             return x[index];
346         } else {
347             return *(last - (n - 1 - index));
348         }
349     }
350
351     /// \endcond
352
353     // state representation: next output is o(x(i))
354     //   x[0]  ... x[k] x[k+1] ... x[n-1]   represents
355     //  x(i-k) ... x(i) x(i+1) ... x(i-k+n-1)
356
357     UIntType x[n]; 
358     std::size_t i;
359 };
360
361 /// \cond show_private
362
363 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
364 //  A definition is required even for integral static constants
365 #define BOOST_RANDOM_MT_DEFINE_CONSTANT(type, name)                         \
366 template<class UIntType, std::size_t w, std::size_t n, std::size_t m,       \
367     std::size_t r, UIntType a, std::size_t u, UIntType d, std::size_t s,    \
368     UIntType b, std::size_t t, UIntType c, std::size_t l, UIntType f>       \
369 const type mersenne_twister_engine<UIntType,w,n,m,r,a,u,d,s,b,t,c,l,f>::name
370 BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, word_size);
371 BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, state_size);
372 BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, shift_size);
373 BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, mask_bits);
374 BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, xor_mask);
375 BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, tempering_u);
376 BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, tempering_d);
377 BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, tempering_s);
378 BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, tempering_b);
379 BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, tempering_t);
380 BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, tempering_c);
381 BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, tempering_l);
382 BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, initialization_multiplier);
383 BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, default_seed);
384 BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, parameter_a);
385 BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, output_u );
386 BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, output_s);
387 BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, output_b);
388 BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, output_t);
389 BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, output_c);
390 BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, output_l);
391 BOOST_RANDOM_MT_DEFINE_CONSTANT(bool, has_fixed_range);
392 #undef BOOST_RANDOM_MT_DEFINE_CONSTANT
393 #endif
394
395 template<class UIntType,
396          std::size_t w, std::size_t n, std::size_t m, std::size_t r,
397          UIntType a, std::size_t u, UIntType d, std::size_t s,
398          UIntType b, std::size_t t,
399          UIntType c, std::size_t l, UIntType f>
400 void
401 mersenne_twister_engine<UIntType,w,n,m,r,a,u,d,s,b,t,c,l,f>::twist()
402 {
403     const UIntType upper_mask = (~static_cast<UIntType>(0)) << r;
404     const UIntType lower_mask = ~upper_mask;
405
406     const std::size_t unroll_factor = 6;
407     const std::size_t unroll_extra1 = (n-m) % unroll_factor;
408     const std::size_t unroll_extra2 = (m-1) % unroll_factor;
409
410     // split loop to avoid costly modulo operations
411     {  // extra scope for MSVC brokenness w.r.t. for scope
412         for(std::size_t j = 0; j < n-m-unroll_extra1; j++) {
413             UIntType y = (x[j] & upper_mask) | (x[j+1] & lower_mask);
414             x[j] = x[j+m] ^ (y >> 1) ^ ((x[j+1]&1) * a);
415         }
416     }
417     {
418         for(std::size_t j = n-m-unroll_extra1; j < n-m; j++) {
419             UIntType y = (x[j] & upper_mask) | (x[j+1] & lower_mask);
420             x[j] = x[j+m] ^ (y >> 1) ^ ((x[j+1]&1) * a);
421         }
422     }
423     {
424         for(std::size_t j = n-m; j < n-1-unroll_extra2; j++) {
425             UIntType y = (x[j] & upper_mask) | (x[j+1] & lower_mask);
426             x[j] = x[j-(n-m)] ^ (y >> 1) ^ ((x[j+1]&1) * a);
427         }
428     }
429     {
430         for(std::size_t j = n-1-unroll_extra2; j < n-1; j++) {
431             UIntType y = (x[j] & upper_mask) | (x[j+1] & lower_mask);
432             x[j] = x[j-(n-m)] ^ (y >> 1) ^ ((x[j+1]&1) * a);
433         }
434     }
435     // last iteration
436     UIntType y = (x[n-1] & upper_mask) | (x[0] & lower_mask);
437     x[n-1] = x[m-1] ^ (y >> 1) ^ ((x[0]&1) * a);
438     i = 0;
439 }
440 /// \endcond
441
442 template<class UIntType,
443          std::size_t w, std::size_t n, std::size_t m, std::size_t r,
444          UIntType a, std::size_t u, UIntType d, std::size_t s,
445          UIntType b, std::size_t t,
446          UIntType c, std::size_t l, UIntType f>
447 inline typename
448 mersenne_twister_engine<UIntType,w,n,m,r,a,u,d,s,b,t,c,l,f>::result_type
449 mersenne_twister_engine<UIntType,w,n,m,r,a,u,d,s,b,t,c,l,f>::operator()()
450 {
451     if(i == n)
452         twist();
453     // Step 4
454     UIntType z = x[i];
455     ++i;
456     z ^= ((z >> u) & d);
457     z ^= ((z << s) & b);
458     z ^= ((z << t) & c);
459     z ^= (z >> l);
460     return z;
461 }
462
463 /**
464  * The specializations \mt11213b and \mt19937 are from
465  *
466  *  @blockquote
467  *  "Mersenne Twister: A 623-dimensionally equidistributed
468  *  uniform pseudo-random number generator", Makoto Matsumoto
469  *  and Takuji Nishimura, ACM Transactions on Modeling and
470  *  Computer Simulation: Special Issue on Uniform Random Number
471  *  Generation, Vol. 8, No. 1, January 1998, pp. 3-30. 
472  *  @endblockquote
473  */
474 typedef mersenne_twister_engine<uint32_t,32,351,175,19,0xccab8ee7,
475     11,0xffffffff,7,0x31b6ab00,15,0xffe50000,17,1812433253> mt11213b;
476
477 /**
478  * The specializations \mt11213b and \mt19937 are from
479  *
480  *  @blockquote
481  *  "Mersenne Twister: A 623-dimensionally equidistributed
482  *  uniform pseudo-random number generator", Makoto Matsumoto
483  *  and Takuji Nishimura, ACM Transactions on Modeling and
484  *  Computer Simulation: Special Issue on Uniform Random Number
485  *  Generation, Vol. 8, No. 1, January 1998, pp. 3-30. 
486  *  @endblockquote
487  */
488 typedef mersenne_twister_engine<uint32_t,32,624,397,31,0x9908b0df,
489     11,0xffffffff,7,0x9d2c5680,15,0xefc60000,18,1812433253> mt19937;
490
491 #if !defined(BOOST_NO_INT64_T) && !defined(BOOST_NO_INTEGRAL_INT64_T)
492 typedef mersenne_twister_engine<uint64_t,64,312,156,31,
493     UINT64_C(0xb5026f5aa96619e9),29,UINT64_C(0x5555555555555555),17,
494     UINT64_C(0x71d67fffeda60000),37,UINT64_C(0xfff7eee000000000),43,
495     UINT64_C(6364136223846793005)> mt19937_64;
496 #endif
497
498 /// \cond show_deprecated
499
500 template<class UIntType,
501          int w, int n, int m, int r,
502          UIntType a, int u, std::size_t s,
503          UIntType b, int t,
504          UIntType c, int l, UIntType v>
505 class mersenne_twister :
506     public mersenne_twister_engine<UIntType,
507         w, n, m, r, a, u, ~(UIntType)0, s, b, t, c, l, 1812433253>
508 {
509     typedef mersenne_twister_engine<UIntType,
510         w, n, m, r, a, u, ~(UIntType)0, s, b, t, c, l, 1812433253> base_type;
511 public:
512     mersenne_twister() {}
513     BOOST_RANDOM_DETAIL_GENERATOR_CONSTRUCTOR(mersenne_twister, Gen, gen)
514     { seed(gen); }
515     BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(mersenne_twister, UIntType, val)
516     { seed(val); }
517     template<class It>
518     mersenne_twister(It& first, It last) : base_type(first, last) {}
519     void seed() { base_type::seed(); }
520     BOOST_RANDOM_DETAIL_GENERATOR_SEED(mersenne_twister, Gen, gen)
521     {
522         detail::generator_seed_seq<Gen> seq(gen);
523         base_type::seed(seq);
524     }
525     BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(mersenne_twister, UIntType, val)
526     { base_type::seed(val); }
527     template<class It>
528     void seed(It& first, It last) { base_type::seed(first, last); }
529 };
530
531 /// \endcond
532
533 } // namespace random
534
535 using random::mt11213b;
536 using random::mt19937;
537 using random::mt19937_64;
538
539 } // namespace boost
540
541 BOOST_RANDOM_PTR_HELPER_SPEC(boost::mt11213b)
542 BOOST_RANDOM_PTR_HELPER_SPEC(boost::mt19937)
543 BOOST_RANDOM_PTR_HELPER_SPEC(boost::mt19937_64)
544
545 #endif // BOOST_RANDOM_MERSENNE_TWISTER_HPP