blitz  Version 1.0.2
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 
67 namespace 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 {
77 public:
78 
79 private:
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.
105  inline twist_int operator() (twist_int s1) {
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 
118 enum { N = 624,
119  PF = 397, // MN98 "m"
120  reference_seed = 4357 };
121 
122  void initialize()
123  {
124  S.resize(N);
125  I = S.end();
126  }
127 
128 public:
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 
148  MersenneTwister(twist_int aa, twist_int bb, twist_int cc) :
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 
161  MersenneTwister(twist_int aa, twist_int bb, twist_int cc,
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.
171  void seed (twist_int seed = reference_seed)
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:
258  State S;
259  State::difference_type I;
260  public:
261  mt_state() { }
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 
292  typedef mt_state T_state;
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 
311 private:
313  const twist_int b_,c_;
314 
315  State S;
316  Iter I;
317 };
318 
319 
323 {
324 public:
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 
331 private:
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
twist_int s0
Definition: mt.h:114
_bz_global blitz::IndexPlaceholder< 9 > r
Definition: indexexpr.h:265
Definition: beta.h:50
twist_int random(void)
Definition: mt.h:241
void reload(void)
Definition: mt.h:218
State::size_type SizeType
Definition: mt.h:82
_bz_global blitz::IndexPlaceholder< 0 > i
Definition: indexexpr.h:256
BitMixer twist_
Definition: mt.h:312
Definition: mt.h:118
State::iterator Iter
Definition: mt.h:83
static const twist_int b_[n]
Definition: mt.h:334
static const twist_int a_[n]
Definition: mt.h:333
_bz_global blitz::IndexPlaceholder< 1 > j
Definition: indexexpr.h:257
MersenneTwister(twist_int aa, twist_int bb, twist_int cc, twist_int initial_seed)
Definition: mt.h:161
MersenneTwister(twist_int aa, twist_int bb, twist_int cc)
Definition: mt.h:148
void seed(State seed_vector)
Definition: mt.h:188
static MersenneTwister create(unsigned int i)
Definition: mt.h:325
twist_int operator()(twist_int s1)
Definition: mt.h:105
This class creates MersenneTwisters with different parameters indexed by and ID number.
Definition: mt.h:322
void setState(const std::string &s)
Definition: mt.h:306
State::difference_type I
Definition: mt.h:259
twist_int low_mask(twist_int s1) const
Definition: mt.h:91
twist_int high_mask(twist_int s1) const
Definition: mt.h:100
unsigned long twist_int
Definition: mt.h:70
static const twist_int c_[n]
Definition: mt.h:335
void setState(const T_state &s)
Definition: mt.h:298
void initialize()
Definition: mt.h:122
mt_state T_state
Definition: mt.h:292
std::string str() const
Definition: mt.h:271
std::string getStateString() const
Definition: mt.h:294
MersenneTwister(twist_int initial_seed)
Definition: mt.h:155
_bz_global blitz::IndexPlaceholder< 2 > k
Definition: indexexpr.h:258
Definition: mt.h:75
std::vector< twist_int > State
Definition: mt.h:81
T_state getState() const
Definition: mt.h:293
BitMixer(twist_int k)
Definition: mt.h:88
#define bool
Definition: compiler.h:100
static const unsigned int n
Definition: mt.h:332
mt_state(const std::string &s)
Definition: mt.h:263
const twist_int c_
Definition: mt.h:313
mt_state()
Definition: mt.h:261
MersenneTwister()
Definition: mt.h:129
_bz_global blitz::IndexPlaceholder< 10 > s
Definition: indexexpr.h:266
Definition: mt.h:119
Iter I
Definition: mt.h:316
BitMixer()
Definition: mt.h:87
void seed(twist_int seed=reference_seed)
Definition: mt.h:171
_bz_global blitz::IndexPlaceholder< 5 > n
Definition: indexexpr.h:261
mt_state(State s, State::difference_type i)
Definition: mt.h:262
const twist_int b_
Definition: mt.h:313
State S
Definition: mt.h:315
const twist_int K
Definition: mt.h:115
State S
Definition: mt.h:258