1 /* boost random/mersenne_twister.hpp header file
2  *
3  * Copyright Jens Maurer 2000-2001
4  * Distributed under the Boost Software License, Version 1.0. (See
5  * accompanying file LICENSE_1_0.txt or copy at
6  * http://www.boost.org/LICENSE_1_0.txt)
7  *
8  * See http://www.boost.org for most recent version including documentation.
9  *
10  * $Id: mersenne_twister.hpp,v 1.20 2005/07/21 22:04:31 jmaurer Exp $
11  *
12  * Revision history
13  *  2001-02-18  moved to individual header files
14 */
15 
16 #ifndef DLIB_BOOST_RANDOM_MERSENNE_TWISTER_HPP
17 #define DLIB_BOOST_RANDOM_MERSENNE_TWISTER_HPP
18 
19 #include <iostream>
20 #include <algorithm>     // std::copy
21 #include <stdexcept>
22 #include "../uintn.h"
23 #include "../serialize.h"
24 
25 namespace dlib
26 {
27     namespace random_helpers
28     {
29 
30     // ------------------------------------------------------------------------------------
31 
32         // http://www.math.keio.ac.jp/matumoto/emt.html
33         template<
34             class UIntType,
35             int w,
36             int n,
37             int m,
38             int r,
39             UIntType a,
40             int u,
41             int s,
42             UIntType b,
43             int t,
44             UIntType c,
45             int l,
46             UIntType val
47             >
48         class mersenne_twister
49         {
50         public:
51             typedef UIntType result_type;
52             const static int word_size = w;
53             const static int state_size = n;
54             const static int shift_size = m;
55             const static int mask_bits = r;
56             const static UIntType parameter_a = a;
57             const static int output_u = u;
58             const static int output_s = s;
59             const static UIntType output_b = b;
60             const static int output_t = t;
61             const static UIntType output_c = c;
62             const static int output_l = l;
63 
64             const static bool has_fixed_range = false;
65 
mersenne_twister()66             mersenne_twister() { seed(); }
67 
mersenne_twister(UIntType value)68             explicit mersenne_twister(UIntType value) { seed(value); }
69 
seed()70             void seed () { seed(UIntType(5489)); }
71 
72             // compiler-generated copy ctor and assignment operator are fine
73 
seed(UIntType value)74             void seed(UIntType value)
75             {
76                 // New seeding algorithm from
77                 // http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html
78                 // In the previous versions, MSBs of the seed affected only MSBs of the
79                 // state x[].
80                 const UIntType mask = ~0u;
81                 x[0] = value & mask;
82                 for (i = 1; i < n; i++) {
83                     // See Knuth "The Art of Computer Programming" Vol. 2, 3rd ed., page 106
84                     x[i] = (1812433253UL * (x[i-1] ^ (x[i-1] >> (w-2))) + i) & mask;
85                 }
86             }
87 
min()88             result_type min() const { return 0; }
max()89             result_type max() const
90             {
91                 // avoid "left shift count >= with of type" warning
92                 result_type res = 0;
93                 for(int i = 0; i < w; ++i)
94                     res |= (1u << i);
95                 return res;
96             }
97 
98             result_type operator()();
99 
serialize(const mersenne_twister & item,std::ostream & out)100             friend void serialize(
101                 const mersenne_twister& item,
102                 std::ostream& out
103             )
104             {
105                 dlib::serialize(item.x, out);
106                 dlib::serialize(item.i, out);
107             }
108 
deserialize(mersenne_twister & item,std::istream & in)109             friend void deserialize(
110                 mersenne_twister& item,
111                 std::istream& in
112             )
113             {
114                 dlib::deserialize(item.x, in);
115                 dlib::deserialize(item.i, in);
116             }
117 
118         private:
119 
120             void twist(int block);
121 
122             // state representation: next output is o(x(i))
123             //   x[0]  ... x[k] x[k+1] ... x[n-1]     x[n]     ... x[2*n-1]   represents
124             //  x(i-k) ... x(i) x(i+1) ... x(i-k+n-1) x(i-k-n) ... x[i(i-k-1)]
125             // The goal is to always have x(i-n) ... x(i-1) available for
126             // operator== and save/restore.
127 
128             UIntType x[2*n];
129             int i;
130         };
131 
132     // ------------------------------------------------------------------------------------
133 
134         template<
135             class UIntType, int w, int n, int m, int r, UIntType a, int u,
136             int s, UIntType b, int t, UIntType c, int l, UIntType val
137             >
twist(int block)138         void mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::twist(
139             int block
140         )
141         {
142             const UIntType upper_mask = (~0u) << r;
143             const UIntType lower_mask = ~upper_mask;
144 
145             if(block == 0) {
146                 for(int j = n; j < 2*n; j++) {
147                     UIntType y = (x[j-n] & upper_mask) | (x[j-(n-1)] & lower_mask);
148                     x[j] = x[j-(n-m)] ^ (y >> 1) ^ (y&1 ? a : 0);
149                 }
150             } else if (block == 1) {
151                 // split loop to avoid costly modulo operations
152                 {  // extra scope for MSVC brokenness w.r.t. for scope
153                     for(int j = 0; j < n-m; j++) {
154                         UIntType y = (x[j+n] & upper_mask) | (x[j+n+1] & lower_mask);
155                         x[j] = x[j+n+m] ^ (y >> 1) ^ (y&1 ? a : 0);
156                     }
157                 }
158 
159                 for(int j = n-m; j < n-1; j++) {
160                     UIntType y = (x[j+n] & upper_mask) | (x[j+n+1] & lower_mask);
161                     x[j] = x[j-(n-m)] ^ (y >> 1) ^ (y&1 ? a : 0);
162                 }
163                 // last iteration
164                 UIntType y = (x[2*n-1] & upper_mask) | (x[0] & lower_mask);
165                 x[n-1] = x[m-1] ^ (y >> 1) ^ (y&1 ? a : 0);
166                 i = 0;
167             }
168         }
169 
170     // ------------------------------------------------------------------------------------
171 
172         template<
173             class UIntType, int w, int n, int m, int r, UIntType a, int u,
174             int s, UIntType b, int t, UIntType c, int l, UIntType val
175             >
176         inline typename mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::result_type
operator()177         mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::operator()(
178         )
179         {
180             if(i == n)
181                 twist(0);
182             else if(i >= 2*n)
183                 twist(1);
184             // Step 4
185             UIntType z = x[i];
186             ++i;
187             z ^= (z >> u);
188             z ^= ((z << s) & b);
189             z ^= ((z << t) & c);
190             z ^= (z >> l);
191             return z;
192         }
193 
194     // ------------------------------------------------------------------------------------
195 
196     } // namespace random
197 
198 
199     typedef random_helpers::mersenne_twister<uint32,32,351,175,19,0xccab8ee7,11,
200     7,0x31b6ab00,15,0xffe50000,17, 0xa37d3c92> mt11213b;
201 
202     // validation by experiment from mt19937.c
203     typedef random_helpers::mersenne_twister<uint32,32,624,397,31,0x9908b0df,11,
204     7,0x9d2c5680,15,0xefc60000,18, 3346425566U> mt19937;
205 
206 } // namespace dlib
207 
208 
209 #endif // DLIB_BOOST_RANDOM_MERSENNE_TWISTER_HPP
210 
211