blitz Version 1.0.2
Loading...
Searching...
No Matches
mt.h
Go to the documentation of this file.
1// -*- C++ -*-
2
3/*
4 * $Id$
5 *
6 * A C-program for MT19937: Integer version (1998/4/6)
7 * genrand() generates one pseudorandom unsigned integer (32bit)
8 * which is uniformly distributed among 0 to 2^32-1 for each
9 * call. sgenrand(seed) set initial values to the working area
10 * of 624 words. Before genrand(), sgenrand(seed) must be
11 * called once. (seed is any 32-bit integer except for 0).
12 * Coded by Takuji Nishimura, considering the suggestions by
13 * Topher Cooper and Marc Rieffel in July-Aug. 1997.
14 *
15 * This library is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU Library General Public
17 * License as published by the Free Software Foundation; either
18 * version 2 of the License, or (at your option) any later
19 * version.
20 * This library is distributed in the hope that it will be useful,
21 * but WITHOUT ANY WARRANTY; without even the implied warranty of
22 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
23 * See the GNU Library General Public License for more details.
24 * You should have received a copy of the GNU Library General
25 * Public License along with this library; if not, write to the
26 * Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
27 * 02111-1307 USA
28 *
29 * Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura.
30 * When you use this, send an email to: matumoto@math.keio.ac.jp
31 * with an appropriate reference to your work.
32 *
33 * REFERENCE
34 * M. Matsumoto and T. Nishimura,
35 * "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform
36 * Pseudo-Random Number Generator",
37 * ACM Transactions on Modeling and Computer Simulation,
38 * Vol. 8, No. 1, January 1998, pp 3--30.
39 *
40 * See
41 * http://www.math.keio.ac.jp/~matumoto/emt.html
42 * and
43 * http://www.acm.org/pubs/citations/journals/tomacs/1998-8-1/p3-matsumoto/
44 *
45 */
46
47#ifndef BZ_RAND_MT
48#define BZ_RAND_MT
49
50#include <blitz/blitz.h>
51
52#include <vector>
53#include <string>
54#include <sstream>
55#include <iostream>
56#include <iterator>
57
58#ifndef UINT_MAX
59 #include <limits.h>
60#endif
61
62#ifdef BZ_HAVE_BOOST_SERIALIZATION
63#include <boost/serialization/serialization.hpp>
64#include <boost/serialization/vector.hpp>
65#endif
66
67namespace ranlib {
68
69#if UINT_MAX < 4294967295U
70 typedef unsigned long twist_int; // must be at least 32 bits
71#else
72 typedef unsigned int twist_int;
73#endif
74
76{
77public:
78
79private:
80
81 typedef std::vector<twist_int> State;
82 typedef State::size_type SizeType;
83 typedef State::iterator Iter;
84
85 // Implements Step 2 and half of Step 3 in the MN98 description
86 struct BitMixer {
87 BitMixer() : s0(0), K(0x9908b0df) {};
88 BitMixer(twist_int k) : s0(0), K(k) {};
89
90 // Return 0 if lsb of s1=0, a if lsb of s1=1.
91 inline twist_int low_mask (twist_int s1) const {
92 // This does not actually result in a branch because it's
93 // equivalent to ( -(s1 & 1u) ) & K
94 return (s1&1u) ? K : 0u;
95 }
96
97 // Return y>>1 in MN98 (step 2 + part of 3).
98 // y = (x[i] AND u) OR (x[i+1 mod n] AND ll), where s0=x[i] and
99 // s1=x[i+1 mod n]
100 inline twist_int high_mask (twist_int s1) const {
101 return ( (s0&0x80000000) | (s1&0x7fffffff) ) >>1;
102 }
103
104 // Calculate (half of) step 3 in MN98.
106 // (y>>1) XOR (0 if lsb of s1=0, a if lsb of s1=1)
107 // x[i+m] is XORed in reload
108 // (Note that it is OK to call low_mask with s1 (x[i+1]) and not
109 // with y, like MN98 does, because s1&1 == y&1 by construction.
110 twist_int r = high_mask(s1) ^ low_mask(s1);
111 s0 = s1;
112 return r;
113 }
114 twist_int s0; // this is "x[i]" in the MN98 description
115 const twist_int K; // MN98 "a" vector
116 };
117
118enum { N = 624,
119 PF = 397, // MN98 "m"
120 reference_seed = 4357 };
121
123 {
124 S.resize(N);
125 I = S.end();
126 }
127
128public:
129 MersenneTwister() : b_(0x9D2C5680), c_(0xEFC60000)
130 {
131 initialize();
132 seed();
133
134 // There is a problem: static initialization + templates do not
135 // mix very well in C++. If you have a static member
136 // of a class template, there is no guarantee on its order iin
137 // static initialization. This MersenneTwister class is used
138 // elsewhere as a static member of a template class, and it is
139 // possible (in fact, I've done so) to create a static initializer
140 // that will invoke the seed() method of this object before its
141 // ctor has been called (result: crash).
142 // ANSI C++ is stranger than fiction.
143 // Currently the documentation forbids using RNGs from
144 // static initializers. There doesn't seem to be a good
145 // fix.
146 }
147
149 twist_(aa), b_(bb), c_(cc)
150 {
151 initialize();
152 seed();
153 }
154
155 MersenneTwister(twist_int initial_seed) : b_(0x9D2C5680), c_(0xEFC60000)
156 {
157 initialize();
158 seed(initial_seed);
159 }
160
162 twist_int initial_seed) : twist_(aa), b_(bb), c_(cc)
163 {
164 initialize();
165 seed(initial_seed);
166 }
167
168 // Seed. Uses updated seed algorithm from mt19937ar. The old
169 // algorithm would yield sequences that were close from seeds that
170 // were close.
172 {
173 // seed cannot equal 0
174 if (seed == 0)
176
177 S[0] = seed & 0xFFFFFFFF;
178 for (SizeType mti=1; mti<S.size(); ++mti) {
179 S[mti] = (1812433253U * (S[mti-1] ^ (S[mti-1] >> 30)) + mti);
180 S[mti] &= 0xffffffffU;
181 }
182
183 reload();
184 }
185
186 // Seed by array, swiped directly from mt19937ar. Gives a larger
187 // initial seed space.
188 void seed (State seed_vector)
189 {
190 SizeType i, j, k;
191 seed(19650218U);
192 i=1; j=0;
193 const SizeType N=S.size();
194 const SizeType n=seed_vector.size();
195 k = (N>n ? N : n);
196 for (; k; k--) {
197 S[i] = (S[i] ^ ((S[i-1] ^ (S[i-1] >> 30)) * 1664525U))
198 + seed_vector[j] + j; /* non linear */
199 S[i] &= 0xffffffffU; /* for WORDSIZE > 32 machines */
200 i++; j++;
201 if (i>=N) { S[0] = S[N-1]; i=1; }
202 if (j>=n) j=0;
203 }
204 for (k=N-1; k; k--) {
205 S[i] = (S[i] ^ ((S[i-1] ^ (S[i-1] >> 30)) * 1566083941UL))
206 - i; /* non linear */
207 S[i] &= 0xffffffffU; /* for WORDSIZE > 32 machines */
208 i++;
209 if (i>=N) { S[0] = S[N-1]; i=1; }
210 }
211
212 S[0] = 0x80000000U; /* MSB is 1; assuring non-zero initial array */
213
214 reload();
215 }
216
217 // generate a full new x array
218 void reload (void)
219 {
220 // This check is required because it is possible to call random()
221 // before the constructor. See the note above about static
222 // initialization.
223
224 Iter p0 = S.begin();
225 Iter pM = p0 + PF;
226 twist_ (S[0]); // set x[i]=x[0] in the twister (prime the pump)
227 for (Iter pf_end = S.begin()+(N-PF); p0 != pf_end; ++p0, ++pM)
228 // mt[kk] = mt[kk+m] XOR ((y>>1)XOR(mag01), as calc by BitMixer)
229 *p0 = *pM ^ twist_ (p0[1]);
230
231 // This is the "modulo part" where kk+m rolls over
232 pM = S.begin();
233 for (Iter s_end = S.begin()+(N-1); p0 != s_end; ++p0, ++pM)
234 *p0 = *pM ^ twist_ (p0[1]);
235 // and final element where kk+1 rolls over
236 *p0 = *pM ^ twist_ (S[0]);
237
238 I = S.begin();
239 }
240
241 inline twist_int random (void)
242 {
243 if (I >= S.end()) reload();
244 // get next word from array
245 twist_int y = *I++;
246 // Step 4+5 in MN98, multiply by tempering matrix
247 y ^= (y >> 11);
248 y ^= (y << 7) & b_;
249 y ^= (y << 15) & c_;
250 y ^= (y >> 18);
251 return y;
252 }
253
254 // functions for getting/setting state
255 class mt_state {
256 friend class MersenneTwister;
257 private:
259 State::difference_type I;
260 public:
262 mt_state(State s, State::difference_type i) : S(s), I(i) { }
263 mt_state(const std::string& s) {
264 std::istringstream is(s);
265 is >> I;
266 S = State(std::istream_iterator<twist_int>(is),
267 std::istream_iterator<twist_int>());
268 assert(!S.empty());
269 }
270 operator bool() const { return !S.empty(); }
271 std::string str() const {
272 if (S.empty())
273 return std::string();
274 std::ostringstream os;
275 os << I << " ";
276 std::copy(S.begin(), S.end(),
277 std::ostream_iterator<twist_int>(os," "));
278 return os.str();
279 }
280#ifdef BZ_HAVE_BOOST_SERIALIZATION
281 friend class boost::serialization::access;
284 template<class T_arch>
285 void serialize(T_arch& ar, const unsigned int version) {
286 ar & S & I;
287 };
288#endif
289
290 };
291
293 T_state getState() const { return T_state(S, I-S.begin()); }
294 std::string getStateString() const {
295 T_state tmp(S, I-S.begin());
296 return tmp.str();
297 }
298 void setState(const T_state& s) {
299 if (!s) {
300 std::cerr << "Error: state is empty" << std::endl;
301 return;
302 }
303 S = s.S;
304 I = S.begin() + s.I;
305 }
306 void setState(const std::string& s) {
307 T_state tmp(s);
308 setState(tmp);
309 }
310
311private:
314
317};
318
319
323{
324public:
325 static MersenneTwister create(unsigned int i) {
326 // We only have n different parameter sets
327 i = i % n;
328 return MersenneTwister(a_[i], b_[i], c_[i]);
329 };
330
331private:
332 static const unsigned int n=48;
333 static const twist_int a_[n];
334 static const twist_int b_[n];
335 static const twist_int c_[n];
336};
337
338}
339
340#endif // BZ_RAND_MT
This class creates MersenneTwisters with different parameters indexed by and ID number.
Definition mt.h:323
static const unsigned int n
Definition mt.h:332
static const twist_int a_[n]
Definition mt.h:333
static MersenneTwister create(unsigned int i)
Definition mt.h:325
static const twist_int c_[n]
Definition mt.h:335
static const twist_int b_[n]
Definition mt.h:334
mt_state(const std::string &s)
Definition mt.h:263
mt_state(State s, State::difference_type i)
Definition mt.h:262
State::difference_type I
Definition mt.h:259
mt_state()
Definition mt.h:261
std::string str() const
Definition mt.h:271
State S
Definition mt.h:258
Definition mt.h:76
MersenneTwister(twist_int initial_seed)
Definition mt.h:155
State S
Definition mt.h:315
void setState(const T_state &s)
Definition mt.h:298
void setState(const std::string &s)
Definition mt.h:306
const twist_int b_
Definition mt.h:313
const twist_int c_
Definition mt.h:313
std::string getStateString() const
Definition mt.h:294
MersenneTwister()
Definition mt.h:129
MersenneTwister(twist_int aa, twist_int bb, twist_int cc)
Definition mt.h:148
void initialize()
Definition mt.h:122
T_state getState() const
Definition mt.h:293
std::vector< twist_int > State
Definition mt.h:81
BitMixer twist_
Definition mt.h:312
mt_state T_state
Definition mt.h:292
void seed(twist_int seed=reference_seed)
Definition mt.h:171
Iter I
Definition mt.h:316
@ PF
Definition mt.h:119
@ N
Definition mt.h:118
@ reference_seed
Definition mt.h:120
State::size_type SizeType
Definition mt.h:82
MersenneTwister(twist_int aa, twist_int bb, twist_int cc, twist_int initial_seed)
Definition mt.h:161
twist_int random(void)
Definition mt.h:241
State::iterator Iter
Definition mt.h:83
void reload(void)
Definition mt.h:218
void seed(State seed_vector)
Definition mt.h:188
#define bool
Definition compiler.h:100
Definition beta.h:50
unsigned long twist_int
Definition mt.h:70
twist_int low_mask(twist_int s1) const
Definition mt.h:91
BitMixer(twist_int k)
Definition mt.h:88
twist_int operator()(twist_int s1)
Definition mt.h:105
twist_int high_mask(twist_int s1) const
Definition mt.h:100
const twist_int K
Definition mt.h:115
twist_int s0
Definition mt.h:114