CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

MTwistEngine.cc
Go to the documentation of this file.
1 // $Id: MTwistEngine.cc,v 1.6 2010/06/16 17:24:53 garren Exp $
2 // -*- C++ -*-
3 //
4 // -----------------------------------------------------------------------
5 // HEP Random
6 // --- MTwistEngine ---
7 // class implementation file
8 // -----------------------------------------------------------------------
9 // A "fast, compact, huge-period generator" based on M. Matsumoto and
10 // T. Nishimura, "Mersenne Twister: A 623-dimensionally equidistributed
11 // uniform pseudorandom number generator", to appear in ACM Trans. on
12 // Modeling and Computer Simulation. It is a twisted GFSR generator
13 // with a Mersenne-prime period of 2^19937-1, uniform on open interval (0,1)
14 // =======================================================================
15 // Ken Smith - Started initial draft: 14th Jul 1998
16 // - Optimized to get std::pow() out of flat() method
17 // - Added conversion operators: 6th Aug 1998
18 // J. Marraffino - Added some explicit casts to deal with
19 // machines where sizeof(int) != sizeof(long) 22 Aug 1998
20 // M. Fischler - Modified constructors such that no two
21 // seeds will match sequences, no single/double
22 // seeds will match, explicit seeds give
23 // determined results, and default will not
24 // match any of the above or other defaults.
25 // - Modified use of the various exponents of 2
26 // to avoid per-instance space overhead and
27 // correct the rounding procedure 16 Sep 1998
28 // J. Marfaffino - Remove dependence on hepString class 13 May 1999
29 // M. Fischler - In restore, checkFile for file not found 03 Dec 2004
30 // M. Fischler - Methods for distrib. instacne save/restore 12/8/04
31 // M. Fischler - split get() into tag validation and
32 // getState() for anonymous restores 12/27/04
33 // M. Fischler - put/get for vectors of ulongs 3/14/05
34 // M. Fischler - State-saving using only ints, for portability 4/12/05
35 // M. Fischler - Improved seeding in setSeed (Savanah bug #17479) 11/15/06
36 // - (Possible optimization - now that the seeding is improved,
37 // is it still necessary for ctor to "warm up" by discarding
38 // 2000 iterations?)
39 //
40 // =======================================================================
41 
42 #include "CLHEP/Random/defs.h"
43 #include "CLHEP/Random/Random.h"
44 #include "CLHEP/Random/MTwistEngine.h"
45 #include "CLHEP/Random/engineIDulong.h"
46 #include <string.h> // for strcmp
47 #include <cstdlib> // for std::abs(int)
48 
49 namespace CLHEP {
50 
51 static const int MarkerLen = 64; // Enough room to hold a begin or end marker.
52 
53 std::string MTwistEngine::name() const {return "MTwistEngine";}
54 
55 int MTwistEngine::numEngines = 0;
56 int MTwistEngine::maxIndex = 215;
57 
60 {
61  int cycle = std::abs(int(numEngines/maxIndex));
62  int curIndex = std::abs(int(numEngines%maxIndex));
63  long mask = ((cycle & 0x007fffff) << 8);
64  long seedlist[2];
65  HepRandom::getTheTableSeeds( seedlist, curIndex );
66  seedlist[0] = (seedlist[0])^mask;
67  seedlist[1] = 0;
68  setSeeds( seedlist, numEngines );
69  count624=0;
70  ++numEngines;
71  for( int i=0; i < 2000; ++i ) flat(); // Warm up just a bit
72 }
73 
76 {
77  long seedlist[2]={seed,17587};
78  setSeeds( seedlist, 0 );
79  count624=0;
80  for( int i=0; i < 2000; ++i ) flat(); // Warm up just a bit
81 }
82 
83 MTwistEngine::MTwistEngine(int rowIndex, int colIndex)
85 {
86  int cycle = std::abs(int(rowIndex/maxIndex));
87  int row = std::abs(int(rowIndex%maxIndex));
88  int col = std::abs(int(colIndex%2));
89  long mask = (( cycle & 0x000007ff ) << 20 );
90  long seedlist[2];
91  HepRandom::getTheTableSeeds( seedlist, row );
92  seedlist[0] = (seedlist[col])^mask;
93  seedlist[1] = 690691;
94  setSeeds(seedlist, 4444772);
95  count624=0;
96  for( int i=0; i < 2000; ++i ) flat(); // Warm up just a bit
97 }
98 
99 MTwistEngine::MTwistEngine( std::istream& is )
100 : HepRandomEngine()
101 {
102  is >> *this;
103 }
104 
106 
108  unsigned int y;
109 
110  if( count624 >= N ) {
111  register int i;
112 
113  for( i=0; i < NminusM; ++i ) {
114  y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
115  mt[i] = mt[i+M] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
116  }
117 
118  for( ; i < N-1 ; ++i ) {
119  y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
120  mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
121  }
122 
123  y = (mt[i] & 0x80000000) | (mt[0] & 0x7fffffff);
124  mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
125 
126  count624 = 0;
127  }
128 
129  y = mt[count624];
130  y ^= ( y >> 11);
131  y ^= ((y << 7 ) & 0x9d2c5680);
132  y ^= ((y << 15) & 0xefc60000);
133  y ^= ( y >> 18);
134 
135  return y * twoToMinus_32() + // Scale to range
136  (mt[count624++] >> 11) * twoToMinus_53() + // fill remaining bits
137  nearlyTwoToMinus_54(); // make sure non-zero
138 }
139 
140 void MTwistEngine::flatArray( const int size, double *vect ) {
141  for( int i=0; i < size; ++i) vect[i] = flat();
142 }
143 
144 void MTwistEngine::setSeed(long seed, int k) {
145 
146  // MF 11/15/06 - Change seeding algorithm to a newer one recommended
147  // by Matsumoto: The original algorithm was
148  // mt[i] = (69069 * mt[i-1]) & 0xffffffff and this gives
149  // problems when the seed bit pattern has lots of zeros
150  // (for example, 0x08000000). Savanah bug #17479.
151 
152  theSeed = seed ? seed : 4357;
153  int mti;
154  const int N1=624;
155  mt[0] = (unsigned int) (theSeed&0xffffffffUL);
156  for (mti=1; mti<N1; mti++) {
157  mt[mti] = (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
158  /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
159  /* In the previous versions, MSBs of the seed affect */
160  /* only MSBs of the array mt[]. */
161  /* 2002/01/09 modified by Makoto Matsumoto */
162  mt[mti] &= 0xffffffffUL;
163  /* for >32 bit machines */
164  }
165  for( int i=1; i < 624; ++i ) {
166  mt[i] ^= k; // MF 9/16/98: distinguish starting points
167  }
168  // MF 11/15/06 This distinction of starting points based on values of k
169  // is kept even though the seeding algorithm itself is improved.
170 }
171 
172 void MTwistEngine::setSeeds(const long *seeds, int k) {
173  setSeed( (*seeds ? *seeds : 43571346), k );
174  int i;
175  for( i=1; i < 624; ++i ) {
176  mt[i] = ( seeds[1] + mt[i] ) & 0xffffffff; // MF 9/16/98
177  }
178  theSeeds = seeds;
179 }
180 
181 void MTwistEngine::saveStatus( const char filename[] ) const
182 {
183  std::ofstream outFile( filename, std::ios::out ) ;
184  if (!outFile.bad()) {
185  outFile << theSeed << std::endl;
186  for (int i=0; i<624; ++i) outFile <<std::setprecision(20) << mt[i] << " ";
187  outFile << std::endl;
188  outFile << count624 << std::endl;
189  }
190 }
191 
192 void MTwistEngine::restoreStatus( const char filename[] )
193 {
194  std::ifstream inFile( filename, std::ios::in);
195  if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) {
196  std::cerr << " -- Engine state remains unchanged\n";
197  return;
198  }
199 
200  if (!inFile.bad() && !inFile.eof()) {
201  inFile >> theSeed;
202  for (int i=0; i<624; ++i) inFile >> mt[i];
203  inFile >> count624;
204  }
205 }
206 
208 {
209  std::cout << std::endl;
210  std::cout << "--------- MTwist engine status ---------" << std::endl;
211  std::cout << std::setprecision(20);
212  std::cout << " Initial seed = " << theSeed << std::endl;
213  std::cout << " Current index = " << count624 << std::endl;
214  std::cout << " Array status mt[] = " << std::endl;
215  for (int i=0; i<624; i+=5) {
216  std::cout << mt[i] << " " << mt[i+1] << " " << mt[i+2] << " "
217  << mt[i+3] << " " << mt[i+4] << std::endl;
218  }
219  std::cout << "----------------------------------------" << std::endl;
220 }
221 
222 MTwistEngine::operator float() {
223  unsigned int y;
224 
225  if( count624 >= N ) {
226  register int i;
227 
228  for( i=0; i < NminusM; ++i ) {
229  y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
230  mt[i] = mt[i+M] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
231  }
232 
233  for( ; i < N-1 ; ++i ) {
234  y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
235  mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
236  }
237 
238  y = (mt[i] & 0x80000000) | (mt[0] & 0x7fffffff);
239  mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
240 
241  count624 = 0;
242  }
243 
244  y = mt[count624++];
245  y ^= ( y >> 11);
246  y ^= ((y << 7 ) & 0x9d2c5680);
247  y ^= ((y << 15) & 0xefc60000);
248  y ^= ( y >> 18);
249 
250  return (float)(y * twoToMinus_32());
251 }
252 
253 MTwistEngine::operator unsigned int() {
254  unsigned int y;
255 
256  if( count624 >= N ) {
257  register int i;
258 
259  for( i=0; i < NminusM; ++i ) {
260  y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
261  mt[i] = mt[i+M] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
262  }
263 
264  for( ; i < N-1 ; ++i ) {
265  y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
266  mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
267  }
268 
269  y = (mt[i] & 0x80000000) | (mt[0] & 0x7fffffff);
270  mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
271 
272  count624 = 0;
273  }
274 
275  y = mt[count624++];
276  y ^= ( y >> 11);
277  y ^= ((y << 7 ) & 0x9d2c5680);
278  y ^= ((y << 15) & 0xefc60000);
279  y ^= ( y >> 18);
280 
281  return y;
282 }
283 
284 std::ostream & MTwistEngine::put ( std::ostream& os ) const
285 {
286  char beginMarker[] = "MTwistEngine-begin";
287  char endMarker[] = "MTwistEngine-end";
288 
289  int pr = os.precision(20);
290  os << " " << beginMarker << " ";
291  os << theSeed << " ";
292  for (int i=0; i<624; ++i) {
293  os << mt[i] << "\n";
294  }
295  os << count624 << " ";
296  os << endMarker << "\n";
297  os.precision(pr);
298  return os;
299 }
300 
301 std::vector<unsigned long> MTwistEngine::put () const {
302  std::vector<unsigned long> v;
303  v.push_back (engineIDulong<MTwistEngine>());
304  for (int i=0; i<624; ++i) {
305  v.push_back(static_cast<unsigned long>(mt[i]));
306  }
307  v.push_back(count624);
308  return v;
309 }
310 
311 std::istream & MTwistEngine::get ( std::istream& is )
312 {
313  char beginMarker [MarkerLen];
314  is >> std::ws;
315  is.width(MarkerLen); // causes the next read to the char* to be <=
316  // that many bytes, INCLUDING A TERMINATION \0
317  // (Stroustrup, section 21.3.2)
318  is >> beginMarker;
319  if (strcmp(beginMarker,"MTwistEngine-begin")) {
320  is.clear(std::ios::badbit | is.rdstate());
321  std::cerr << "\nInput stream mispositioned or"
322  << "\nMTwistEngine state description missing or"
323  << "\nwrong engine type found." << std::endl;
324  return is;
325  }
326  return getState(is);
327 }
328 
329 std::string MTwistEngine::beginTag ( ) {
330  return "MTwistEngine-begin";
331 }
332 
333 std::istream & MTwistEngine::getState ( std::istream& is )
334 {
335  char endMarker [MarkerLen];
336  is >> theSeed;
337  for (int i=0; i<624; ++i) is >> mt[i];
338  is >> count624;
339  is >> std::ws;
340  is.width(MarkerLen);
341  is >> endMarker;
342  if (strcmp(endMarker,"MTwistEngine-end")) {
343  is.clear(std::ios::badbit | is.rdstate());
344  std::cerr << "\nMTwistEngine state description incomplete."
345  << "\nInput stream is probably mispositioned now." << std::endl;
346  return is;
347  }
348  return is;
349 }
350 
351 bool MTwistEngine::get (const std::vector<unsigned long> & v) {
352  if ((v[0] & 0xffffffffUL) != engineIDulong<MTwistEngine>()) {
353  std::cerr <<
354  "\nMTwistEngine get:state vector has wrong ID word - state unchanged\n";
355  return false;
356  }
357  return getState(v);
358 }
359 
360 bool MTwistEngine::getState (const std::vector<unsigned long> & v) {
361  if (v.size() != VECTOR_STATE_SIZE ) {
362  std::cerr <<
363  "\nMTwistEngine get:state vector has wrong length - state unchanged\n";
364  return false;
365  }
366  for (int i=0; i<624; ++i) {
367  mt[i]=v[i+1];
368  }
369  count624 = v[625];
370  return true;
371 }
372 
373 } // namespace CLHEP