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