blitz  Version 0.9
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
rand-mt.h
Go to the documentation of this file.
1 /* A C-program for MT19937: Integer version (1998/4/6) */
2 /* genrand() generates one pseudorandom unsigned integer (32bit) */
3 /* which is uniformly distributed among 0 to 2^32-1 for each */
4 /* call. sgenrand(seed) set initial values to the working area */
5 /* of 624 words. Before genrand(), sgenrand(seed) must be */
6 /* called once. (seed is any 32-bit integer except for 0). */
7 /* Coded by Takuji Nishimura, considering the suggestions by */
8 /* Topher Cooper and Marc Rieffel in July-Aug. 1997. */
9 
10 /* This library is free software; you can redistribute it and/or */
11 /* modify it under the terms of the GNU Library General Public */
12 /* License as published by the Free Software Foundation; either */
13 /* version 2 of the License, or (at your option) any later */
14 /* version. */
15 /* This library is distributed in the hope that it will be useful, */
16 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
17 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. */
18 /* See the GNU Library General Public License for more details. */
19 /* You should have received a copy of the GNU Library General */
20 /* Public License along with this library; if not, write to the */
21 /* Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA */
22 /* 02111-1307 USA */
23 
24 /* Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura. */
25 /* When you use this, send an email to: matumoto@math.keio.ac.jp */
26 /* with an appropriate reference to your work. */
27 
28 /* REFERENCE */
29 /* M. Matsumoto and T. Nishimura, */
30 /* "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform */
31 /* Pseudo-Random Number Generator", */
32 /* ACM Transactions on Modeling and Computer Simulation, */
33 /* Vol. 8, No. 1, January 1998, pp 3--30. */
34 
35 // See http://www.math.keio.ac.jp/~matumoto/emt.html
36 
37 // 1999-01-25 adapted to STL-like idiom
38 // allan@stokes.ca (Allan Stokes) www.stokes.ca
39 
40 #ifndef BZ_RAND_MT
41 #define BZ_RAND_MT
42 
43 #ifndef BZ_BLITZ_H
44  #include <blitz/blitz.h>
45 #endif
46 
47 #include <vector>
48 
49 BZ_NAMESPACE(blitz)
50 
51 // decomposition issues:
52 // machine representation of integer types
53 // output buffer option verses inline post-conditioning
54 
56 {
57 private:
58  typedef unsigned int twist_int; // must be at least 32 bits
59  // larger might be faster
60  typedef vector<twist_int> State;
61  typedef State::iterator Iter;
62 
63  struct BitMixer {
64  enum { K = 0x9908b0df };
65  BitMixer() : s0(0) {}
66  inline friend twist_int low_mask (twist_int s1) {
67  return (s1&1u) ? K : 0u;
68  }
69  inline twist_int high_mask (twist_int s1) const {
70  return ((s0&0x80000000)|(s1&0x7fffffff))>>1;
71  }
72  inline twist_int operator() (twist_int s1) {
73  twist_int r = high_mask(s1) ^ low_mask(s1);
74  s0 = s1;
75  return r;
76  }
77  twist_int s0;
78  };
79 
80 enum { N = 624, PF = 397, reference_seed = 4357 };
81 
82 public:
83  MersenneTwister () {} // S empty will trigger auto-seed
84 
85  void seed (twist_int seed = reference_seed)
86  {
87  if (!S.size()) S.resize(N);
88  enum { Knuth_A = 69069 };
89  twist_int x = seed & 0xFFFFFFFF;
90  Iter s = S.begin();
91  twist_int mask = (seed == reference_seed) ? 0 : 0xFFFFFFFF;
92  for (int j = 0; j < N; ++j) {
93  // adding j here avoids the risk of all zeros
94  // we suppress this term in "compatibility" mode
95  *s++ = (x + (mask & j)) & 0xFFFFFFFF;
96  x *= Knuth_A;
97  }
98  reload();
99  }
100 
101  void reload (void)
102  {
103  if (!S.size()) seed (); // auto-seed detection
104 
105  Iter p0 = S.begin();
106  Iter pM = p0 + PF;
107  BitMixer twist;
108  twist (S[0]); // prime the pump
109  for (Iter pf_end = S.begin()+(N-PF); p0 != pf_end; ++p0, ++pM)
110  *p0 = *pM ^ twist (p0[1]);
111  pM = S.begin();
112  for (Iter s_end = S.begin()+(N-1); p0 != s_end; ++p0, ++pM)
113  *p0 = *pM ^ twist (p0[1]);
114  *p0 = *pM ^ twist (S[0]);
115 
116  I = S.begin();
117  }
118 
119  inline twist_int random (void)
120  {
121  if (I >= S.end()) reload();
122  twist_int y = *I++;
123  y ^= (y >> 11);
124  y ^= (y << 7) & 0x9D2C5680;
125  y ^= (y << 15) & 0xEFC60000;
126  y ^= (y >> 18);
127  return y;
128  }
129 
130 private:
131  State S;
132  Iter I;
133 };
134 
135 
136 // This version returns a double in the range [0,1).
137 
139 
140 public:
142  {
143  // f = 1/(2^32);
144  f = (1.0 / 65536) / 65536;
145  }
146 
147  void randomize(unsigned int seed)
148  {
149  gen_.seed(seed);
150  }
151 
152  double random()
153  {
154  unsigned long y1 = gen_.random();
155  unsigned long y2 = gen_.random();
156 
157  return ((y1 * f) * y2 * f);
158  }
159 
160 private:
162  double f;
163 };
164 
166 
167 #endif // BZ_RAND_MT