blitz  Version 0.9
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
mt.h
Go to the documentation of this file.
1 
2 /*
3  * $Id: mt.h,v 1.6 2005/10/13 20:08:53 julianc Exp $
4  *
5  * A C-program for MT19937: Integer version (1998/4/6)
6  * genrand() generates one pseudorandom unsigned integer (32bit)
7  * which is uniformly distributed among 0 to 2^32-1 for each
8  * call. sgenrand(seed) set initial values to the working area
9  * of 624 words. Before genrand(), sgenrand(seed) must be
10  * called once. (seed is any 32-bit integer except for 0).
11  * Coded by Takuji Nishimura, considering the suggestions by
12  * Topher Cooper and Marc Rieffel in July-Aug. 1997.
13  *
14  * This library is free software; you can redistribute it and/or
15  * modify it under the terms of the GNU Library General Public
16  * License as published by the Free Software Foundation; either
17  * version 2 of the License, or (at your option) any later
18  * version.
19  * This library is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
22  * See the GNU Library General Public License for more details.
23  * You should have received a copy of the GNU Library General
24  * Public License along with this library; if not, write to the
25  * Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
26  * 02111-1307 USA
27  *
28  * Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura.
29  * When you use this, send an email to: matumoto@math.keio.ac.jp
30  * with an appropriate reference to your work.
31  *
32  * REFERENCE
33  * M. Matsumoto and T. Nishimura,
34  * "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform
35  * Pseudo-Random Number Generator",
36  * ACM Transactions on Modeling and Computer Simulation,
37  * Vol. 8, No. 1, January 1998, pp 3--30.
38  *
39  * See
40  * http://www.math.keio.ac.jp/~matumoto/emt.html
41  * and
42  * http://www.acm.org/pubs/citations/journals/tomacs/1998-8-1/p3-matsumoto/
43  *
44  */
45 
46 #ifndef BZ_RAND_MT
47 #define BZ_RAND_MT
48 
49 #include <blitz/blitz.h>
50 
51 #include <vector>
52 #include <string>
53 #include <sstream>
54 #include <iostream>
55 #include <iterator>
56 
57 #ifndef UINT_MAX
58  #include <limits.h>
59 #endif
60 
61 BZ_NAMESPACE(ranlib)
62 
63 class MersenneTwister
64 {
65 public:
66 
67 #if UINT_MAX < 4294967295U
68  typedef unsigned long twist_int; // must be at least 32 bits
69 #else
70  typedef unsigned int twist_int;
71 #endif
72 
73 private:
74 
75 #if defined(BZ_HAVE_NAMESPACES) && defined(BZ_HAVE_STD)
76  typedef std::vector<twist_int> State;
77 #else
78  typedef vector<twist_int> State;
79 #endif
80 
81  typedef State::iterator Iter;
82 
83  struct BitMixer {
84  enum { K = 0x9908b0df };
85  BitMixer() : s0(0) {}
86  inline twist_int low_mask (twist_int s1) const {
87  return (s1&1u) ? K : 0u;
88  }
89  inline twist_int high_mask (twist_int s1) const {
90  return ((s0&0x80000000)|(s1&0x7fffffff))>>1;
91  }
92  inline twist_int operator() (twist_int s1) {
93  twist_int r = high_mask(s1) ^ low_mask(s1);
94  s0 = s1;
95  return r;
96  }
97  twist_int s0;
98  };
99 
100 enum { N = 624, PF = 397, reference_seed = 4357 };
101 
102  void initialize()
103  {
104  S.resize(N);
105  I = S.end();
106  }
107 
108 public:
109  MersenneTwister()
110  {
111  initialize();
112  seed();
113 
114  // There is a problem: static initialization + templates do not
115  // mix very well in C++. If you have a static member
116  // of a class template, there is no guarantee on its order iin
117  // static initialization. This MersenneTwister class is used
118  // elsewhere as a static member of a template class, and it is
119  // possible (in fact, I've done so) to create a static initializer
120  // that will invoke the seed() method of this object before its
121  // ctor has been called (result: crash).
122  // ANSI C++ is stranger than fiction.
123  // Currently the documentation forbids using RNGs from
124  // static initializers. There doesn't seem to be a good
125  // fix.
126  }
127 
128  MersenneTwister(twist_int initial_seed)
129  {
130  initialize();
131  seed(initial_seed);
132  }
133 
134  void seed (twist_int seed = reference_seed)
135  {
136  // seed cannot equal 0
137  if (seed == 0)
138  seed = reference_seed;
139 
140  enum { Knuth_A = 69069 };
141  twist_int x = seed & 0xFFFFFFFF;
142  Iter s = S.begin();
143  twist_int mask = (seed == reference_seed) ? 0 : 0xFFFFFFFF;
144  for (int j = 0; j < N; ++j) {
145  // adding j here avoids the risk of all zeros
146  // we suppress this term in "compatibility" mode
147  *s++ = (x + (mask & j)) & 0xFFFFFFFF;
148  x *= Knuth_A;
149  }
150 
151  reload();
152  }
153 
154  void reload (void)
155  {
156  // This check is required because it is possible to call random()
157  // before the constructor. See the note above about static
158  // initialization.
159 
160  Iter p0 = S.begin();
161  Iter pM = p0 + PF;
162  BitMixer twist;
163  twist (S[0]); // prime the pump
164  for (Iter pf_end = S.begin()+(N-PF); p0 != pf_end; ++p0, ++pM)
165  *p0 = *pM ^ twist (p0[1]);
166  pM = S.begin();
167  for (Iter s_end = S.begin()+(N-1); p0 != s_end; ++p0, ++pM)
168  *p0 = *pM ^ twist (p0[1]);
169  *p0 = *pM ^ twist (S[0]);
170 
171  I = S.begin();
172  }
173 
174  inline twist_int random (void)
175  {
176  if (I >= S.end()) reload();
177  twist_int y = *I++;
178  y ^= (y >> 11);
179  y ^= (y << 7) & 0x9D2C5680;
180  y ^= (y << 15) & 0xEFC60000;
181  y ^= (y >> 18);
182  return y;
183  }
184 
185  // functions for getting/setting state
186  class mt_state {
187  friend class MersenneTwister;
188  private:
189  State S;
190  int I;
191  public:
192  mt_state() { }
193  mt_state(State s, int i) : S(s), I(i) { }
194  mt_state(const std::string& s) {
195  std::istringstream is(s);
196  is >> I;
197  S = State(std::istream_iterator<twist_int>(is),
198  std::istream_iterator<twist_int>());
199  assert(!S.empty());
200  }
201  operator bool() const { return !S.empty(); }
202  std::string str() const {
203  if (S.empty())
204  return std::string();
205  std::ostringstream os;
206  os << I << " ";
207  std::copy(S.begin(), S.end(),
208  std::ostream_iterator<twist_int>(os," "));
209  return os.str();
210  }
211  };
212 
213  typedef mt_state T_state;
214  T_state getState() const { return T_state(S, I-S.begin()); }
215  std::string getStateString() const {
216  T_state tmp(S, I-S.begin());
217  return tmp.str();
218  }
219  void setState(const T_state& s) {
220  if (!s) {
221  std::cerr << "Error: state is empty" << std::endl;
222  return;
223  }
224  S = s.S;
225  I = S.begin() + s.I;
226  }
227  void setState(const std::string& s) {
228  T_state tmp(s);
229  setState(tmp);
230  }
231 
232 private:
233  State S;
234  Iter I;
235 };
236 
237 
239 
240 #endif // BZ_RAND_MT