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