Geant4  10.01.p02
MTwistEngine.cc
Go to the documentation of this file.
1 // $Id:$
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/Random.h"
43 #include "CLHEP/Random/MTwistEngine.h"
44 #include "CLHEP/Random/engineIDulong.h"
45 #include "CLHEP/Utility/atomic_int.h"
46 
47 #include <string.h> // for strcmp
48 #include <cstdlib> // for std::abs(int)
49 
50 namespace CLHEP {
51 
52 namespace {
53  // Number of instances with automatic seed selection
54  CLHEP_ATOMIC_INT_TYPE numberOfEngines(0);
55 
56  // Maximum index into the seed table
57  const int maxIndex = 215;
58 }
59 
60 static const int MarkerLen = 64; // Enough room to hold a begin or end marker.
61 
62 std::string MTwistEngine::name() const {return "MTwistEngine";}
63 
64 MTwistEngine::MTwistEngine()
65 : HepRandomEngine()
66 {
67  int numEngines = numberOfEngines++;
68  int cycle = std::abs(int(numEngines/maxIndex));
69  int curIndex = std::abs(int(numEngines%maxIndex));
70  long mask = ((cycle & 0x007fffff) << 8);
71  long seedlist[2];
72  HepRandom::getTheTableSeeds( seedlist, curIndex );
73  seedlist[0] = (seedlist[0])^mask;
74  seedlist[1] = 0;
75  setSeeds( seedlist, numEngines );
76  count624=0;
77 
78  for( int i=0; i < 2000; ++i ) flat(); // Warm up just a bit
79 }
80 
81 MTwistEngine::MTwistEngine(long seed)
82 : HepRandomEngine()
83 {
84  long seedlist[2]={seed,17587};
85  setSeeds( seedlist, 0 );
86  count624=0;
87  for( int i=0; i < 2000; ++i ) flat(); // Warm up just a bit
88 }
89 
90 MTwistEngine::MTwistEngine(int rowIndex, int colIndex)
91 : HepRandomEngine()
92 {
93  int cycle = std::abs(int(rowIndex/maxIndex));
94  int row = std::abs(int(rowIndex%maxIndex));
95  int col = std::abs(int(colIndex%2));
96  long mask = (( cycle & 0x000007ff ) << 20 );
97  long seedlist[2];
98  HepRandom::getTheTableSeeds( seedlist, row );
99  seedlist[0] = (seedlist[col])^mask;
100  seedlist[1] = 690691;
101  setSeeds(seedlist, 4444772);
102  count624=0;
103  for( int i=0; i < 2000; ++i ) flat(); // Warm up just a bit
104 }
105 
106 MTwistEngine::MTwistEngine( std::istream& is )
107 : HepRandomEngine()
108 {
109  is >> *this;
110 }
111 
112 MTwistEngine::~MTwistEngine() {}
113 
114 double MTwistEngine::flat() {
115  unsigned int y;
116 
117  if( count624 >= N ) {
118  int i;
119 
120  for( i=0; i < NminusM; ++i ) {
121  y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
122  mt[i] = mt[i+M] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
123  }
124 
125  for( ; i < N-1 ; ++i ) {
126  y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
127  mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
128  }
129 
130  y = (mt[i] & 0x80000000) | (mt[0] & 0x7fffffff);
131  mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
132 
133  count624 = 0;
134  }
135 
136  y = mt[count624];
137  y ^= ( y >> 11);
138  y ^= ((y << 7 ) & 0x9d2c5680);
139  y ^= ((y << 15) & 0xefc60000);
140  y ^= ( y >> 18);
141 
142  return y * twoToMinus_32() + // Scale to range
143  (mt[count624++] >> 11) * twoToMinus_53() + // fill remaining bits
144  nearlyTwoToMinus_54(); // make sure non-zero
145 }
146 
147 void MTwistEngine::flatArray( const int size, double *vect ) {
148  for( int i=0; i < size; ++i) vect[i] = flat();
149 }
150 
151 void MTwistEngine::setSeed(long seed, int k) {
152 
153  // MF 11/15/06 - Change seeding algorithm to a newer one recommended
154  // by Matsumoto: The original algorithm was
155  // mt[i] = (69069 * mt[i-1]) & 0xffffffff and this gives
156  // problems when the seed bit pattern has lots of zeros
157  // (for example, 0x08000000). Savanah bug #17479.
158 
159  theSeed = seed ? seed : 4357;
160  int mti;
161  const int N1=624;
162  mt[0] = (unsigned int) (theSeed&0xffffffffUL);
163  for (mti=1; mti<N1; mti++) {
164  mt[mti] = (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
165  /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
166  /* In the previous versions, MSBs of the seed affect */
167  /* only MSBs of the array mt[]. */
168  /* 2002/01/09 modified by Makoto Matsumoto */
169  mt[mti] &= 0xffffffffUL;
170  /* for >32 bit machines */
171  }
172  for( int i=1; i < 624; ++i ) {
173  mt[i] ^= k; // MF 9/16/98: distinguish starting points
174  }
175  // MF 11/15/06 This distinction of starting points based on values of k
176  // is kept even though the seeding algorithm itself is improved.
177 }
178 
179 void MTwistEngine::setSeeds(const long *seeds, int k) {
180  setSeed( (*seeds ? *seeds : 43571346), k );
181  int i;
182  for( i=1; i < 624; ++i ) {
183  mt[i] = ( seeds[1] + mt[i] ) & 0xffffffff; // MF 9/16/98
184  }
185  theSeeds = seeds;
186 }
187 
188 void MTwistEngine::saveStatus( const char filename[] ) const
189 {
190  std::ofstream outFile( filename, std::ios::out ) ;
191  if (!outFile.bad()) {
192  outFile << theSeed << std::endl;
193  for (int i=0; i<624; ++i) outFile <<std::setprecision(20) << mt[i] << " ";
194  outFile << std::endl;
195  outFile << count624 << std::endl;
196  }
197 }
198 
199 void MTwistEngine::restoreStatus( const char filename[] )
200 {
201  std::ifstream inFile( filename, std::ios::in);
202  if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) {
203  std::cerr << " -- Engine state remains unchanged\n";
204  return;
205  }
206 
207  if (!inFile.bad() && !inFile.eof()) {
208  inFile >> theSeed;
209  for (int i=0; i<624; ++i) inFile >> mt[i];
210  inFile >> count624;
211  }
212 }
213 
214 void MTwistEngine::showStatus() const
215 {
216  std::cout << std::endl;
217  std::cout << "--------- MTwist engine status ---------" << std::endl;
218  std::cout << std::setprecision(20);
219  std::cout << " Initial seed = " << theSeed << std::endl;
220  std::cout << " Current index = " << count624 << std::endl;
221  std::cout << " Array status mt[] = " << std::endl;
222  // 2014/06/06 L Garren
223  // the final line has 4 elements, not 5
224  for (int i=0; i<620; i+=5) {
225  std::cout << mt[i] << " " << mt[i+1] << " " << mt[i+2] << " "
226  << mt[i+3] << " " << mt[i+4] << "\n";
227  }
228  std::cout << mt[620] << " " << mt[621] << " " << mt[622] << " "
229  << mt[623] << std::endl;
230  std::cout << "----------------------------------------" << std::endl;
231 }
232 
233 MTwistEngine::operator float() {
234  unsigned int y;
235 
236  if( count624 >= N ) {
237  int i;
238 
239  for( i=0; i < NminusM; ++i ) {
240  y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
241  mt[i] = mt[i+M] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
242  }
243 
244  for( ; i < N-1 ; ++i ) {
245  y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
246  mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
247  }
248 
249  y = (mt[i] & 0x80000000) | (mt[0] & 0x7fffffff);
250  mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
251 
252  count624 = 0;
253  }
254 
255  y = mt[count624++];
256  y ^= ( y >> 11);
257  y ^= ((y << 7 ) & 0x9d2c5680);
258  y ^= ((y << 15) & 0xefc60000);
259  y ^= ( y >> 18);
260 
261  return (float)(y * twoToMinus_32());
262 }
263 
264 MTwistEngine::operator unsigned int() {
265  unsigned int y;
266 
267  if( count624 >= N ) {
268  int i;
269 
270  for( i=0; i < NminusM; ++i ) {
271  y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
272  mt[i] = mt[i+M] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
273  }
274 
275  for( ; i < N-1 ; ++i ) {
276  y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
277  mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
278  }
279 
280  y = (mt[i] & 0x80000000) | (mt[0] & 0x7fffffff);
281  mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
282 
283  count624 = 0;
284  }
285 
286  y = mt[count624++];
287  y ^= ( y >> 11);
288  y ^= ((y << 7 ) & 0x9d2c5680);
289  y ^= ((y << 15) & 0xefc60000);
290  y ^= ( y >> 18);
291 
292  return y;
293 }
294 
295 std::ostream & MTwistEngine::put ( std::ostream& os ) const
296 {
297  char beginMarker[] = "MTwistEngine-begin";
298  char endMarker[] = "MTwistEngine-end";
299 
300  int pr = os.precision(20);
301  os << " " << beginMarker << " ";
302  os << theSeed << " ";
303  for (int i=0; i<624; ++i) {
304  os << mt[i] << "\n";
305  }
306  os << count624 << " ";
307  os << endMarker << "\n";
308  os.precision(pr);
309  return os;
310 }
311 
312 std::vector<unsigned long> MTwistEngine::put () const {
313  std::vector<unsigned long> v;
314  v.push_back (engineIDulong<MTwistEngine>());
315  for (int i=0; i<624; ++i) {
316  v.push_back(static_cast<unsigned long>(mt[i]));
317  }
318  v.push_back(count624);
319  return v;
320 }
321 
322 std::istream & MTwistEngine::get ( std::istream& is )
323 {
324  char beginMarker [MarkerLen];
325  is >> std::ws;
326  is.width(MarkerLen); // causes the next read to the char* to be <=
327  // that many bytes, INCLUDING A TERMINATION \0
328  // (Stroustrup, section 21.3.2)
329  is >> beginMarker;
330  if (strcmp(beginMarker,"MTwistEngine-begin")) {
331  is.clear(std::ios::badbit | is.rdstate());
332  std::cerr << "\nInput stream mispositioned or"
333  << "\nMTwistEngine state description missing or"
334  << "\nwrong engine type found." << std::endl;
335  return is;
336  }
337  return getState(is);
338 }
339 
340 std::string MTwistEngine::beginTag ( ) {
341  return "MTwistEngine-begin";
342 }
343 
344 std::istream & MTwistEngine::getState ( std::istream& is )
345 {
346  char endMarker [MarkerLen];
347  is >> theSeed;
348  for (int i=0; i<624; ++i) is >> mt[i];
349  is >> count624;
350  is >> std::ws;
351  is.width(MarkerLen);
352  is >> endMarker;
353  if (strcmp(endMarker,"MTwistEngine-end")) {
354  is.clear(std::ios::badbit | is.rdstate());
355  std::cerr << "\nMTwistEngine state description incomplete."
356  << "\nInput stream is probably mispositioned now." << std::endl;
357  return is;
358  }
359  return is;
360 }
361 
362 bool MTwistEngine::get (const std::vector<unsigned long> & v) {
363  if ((v[0] & 0xffffffffUL) != engineIDulong<MTwistEngine>()) {
364  std::cerr <<
365  "\nMTwistEngine get:state vector has wrong ID word - state unchanged\n";
366  return false;
367  }
368  return getState(v);
369 }
370 
371 bool MTwistEngine::getState (const std::vector<unsigned long> & v) {
372  if (v.size() != VECTOR_STATE_SIZE ) {
373  std::cerr <<
374  "\nMTwistEngine get:state vector has wrong length - state unchanged\n";
375  return false;
376  }
377  for (int i=0; i<624; ++i) {
378  mt[i]=v[i+1];
379  }
380  count624 = v[625];
381  return true;
382 }
383 
384 } // namespace CLHEP
static const int MarkerLen
Definition: DualRand.cc:67
static ush mask[]
Definition: csz_inflate.cc:332
G4String name
Definition: TRTMaterials.hh:40
double flat()
Definition: G4AblaRandom.cc:47
void setSeeds(const SeedVector &sv)
Set the seeds of the current generator.
Definition: G4INCLRandom.cc:85