Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
MTwistEngine.cc
이 파일의 문서화 페이지로 가기
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"
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 
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 
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)
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 
113 
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 
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 double() {
234  return flat();
235 }
236 
237 MTwistEngine::operator float() {
238  unsigned int y;
239 
240  if( count624 >= N ) {
241  int i;
242 
243  for( i=0; i < NminusM; ++i ) {
244  y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
245  mt[i] = mt[i+M] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
246  }
247 
248  for( ; i < N-1 ; ++i ) {
249  y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
250  mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
251  }
252 
253  y = (mt[i] & 0x80000000) | (mt[0] & 0x7fffffff);
254  mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
255 
256  count624 = 0;
257  }
258 
259  y = mt[count624++];
260  y ^= ( y >> 11);
261  y ^= ((y << 7 ) & 0x9d2c5680);
262  y ^= ((y << 15) & 0xefc60000);
263  y ^= ( y >> 18);
264 
265  return (float)(y * twoToMinus_32());
266 }
267 
268 MTwistEngine::operator unsigned int() {
269  unsigned int y;
270 
271  if( count624 >= N ) {
272  int i;
273 
274  for( i=0; i < NminusM; ++i ) {
275  y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
276  mt[i] = mt[i+M] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
277  }
278 
279  for( ; i < N-1 ; ++i ) {
280  y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
281  mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
282  }
283 
284  y = (mt[i] & 0x80000000) | (mt[0] & 0x7fffffff);
285  mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
286 
287  count624 = 0;
288  }
289 
290  y = mt[count624++];
291  y ^= ( y >> 11);
292  y ^= ((y << 7 ) & 0x9d2c5680);
293  y ^= ((y << 15) & 0xefc60000);
294  y ^= ( y >> 18);
295 
296  return y;
297 }
298 
299 std::ostream & MTwistEngine::put ( std::ostream& os ) const
300 {
301  char beginMarker[] = "MTwistEngine-begin";
302  char endMarker[] = "MTwistEngine-end";
303 
304  int pr = os.precision(20);
305  os << " " << beginMarker << " ";
306  os << theSeed << " ";
307  for (int i=0; i<624; ++i) {
308  os << mt[i] << "\n";
309  }
310  os << count624 << " ";
311  os << endMarker << "\n";
312  os.precision(pr);
313  return os;
314 }
315 
316 std::vector<unsigned long> MTwistEngine::put () const {
317  std::vector<unsigned long> v;
318  v.push_back (engineIDulong<MTwistEngine>());
319  for (int i=0; i<624; ++i) {
320  v.push_back(static_cast<unsigned long>(mt[i]));
321  }
322  v.push_back(count624);
323  return v;
324 }
325 
326 std::istream & MTwistEngine::get ( std::istream& is )
327 {
328  char beginMarker [MarkerLen];
329  is >> std::ws;
330  is.width(MarkerLen); // causes the next read to the char* to be <=
331  // that many bytes, INCLUDING A TERMINATION \0
332  // (Stroustrup, section 21.3.2)
333  is >> beginMarker;
334  if (strcmp(beginMarker,"MTwistEngine-begin")) {
335  is.clear(std::ios::badbit | is.rdstate());
336  std::cerr << "\nInput stream mispositioned or"
337  << "\nMTwistEngine state description missing or"
338  << "\nwrong engine type found." << std::endl;
339  return is;
340  }
341  return getState(is);
342 }
343 
344 std::string MTwistEngine::beginTag ( ) {
345  return "MTwistEngine-begin";
346 }
347 
348 std::istream & MTwistEngine::getState ( std::istream& is )
349 {
350  char endMarker [MarkerLen];
351  is >> theSeed;
352  for (int i=0; i<624; ++i) is >> mt[i];
353  is >> count624;
354  is >> std::ws;
355  is.width(MarkerLen);
356  is >> endMarker;
357  if (strcmp(endMarker,"MTwistEngine-end")) {
358  is.clear(std::ios::badbit | is.rdstate());
359  std::cerr << "\nMTwistEngine state description incomplete."
360  << "\nInput stream is probably mispositioned now." << std::endl;
361  return is;
362  }
363  return is;
364 }
365 
366 bool MTwistEngine::get (const std::vector<unsigned long> & v) {
367  if ((v[0] & 0xffffffffUL) != engineIDulong<MTwistEngine>()) {
368  std::cerr <<
369  "\nMTwistEngine get:state vector has wrong ID word - state unchanged\n";
370  return false;
371  }
372  return getState(v);
373 }
374 
375 bool MTwistEngine::getState (const std::vector<unsigned long> & v) {
376  if (v.size() != VECTOR_STATE_SIZE ) {
377  std::cerr <<
378  "\nMTwistEngine get:state vector has wrong length - state unchanged\n";
379  return false;
380  }
381  for (int i=0; i<624; ++i) {
382  mt[i]=v[i+1];
383  }
384  count624 = v[625];
385  return true;
386 }
387 
388 } // namespace CLHEP
void showStatus() const
static ush mask[]
Definition: csz_inflate.cc:317
Float_t x1[n_points_granero]
Definition: compare.C:5
Float_t y
Definition: compare.C:6
static double nearlyTwoToMinus_54()
static std::string engineName()
Definition: MTwistEngine.h:78
void setSeed(long seed, int)
static double twoToMinus_32()
static const int MarkerLen
Definition: DualRand.cc:67
**D E S C R I P T I O N
virtual std::istream & getState(std::istream &is)
long seed
Definition: chem4.cc:68
Int_t col[ntarg]
Definition: Style.C:29
#define CLHEP_ATOMIC_INT_TYPE
Definition: atomic_int.h:20
void saveStatus(const char filename[]="MTwist.conf") const
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
void restoreStatus(const char filename[]="MTwist.conf")
unsigned int mt[624]
Definition: MTwistEngine.h:88
static const unsigned int VECTOR_STATE_SIZE
Definition: MTwistEngine.h:84
virtual std::istream & get(std::istream &is)
std::string name() const
Definition: MTwistEngine.cc:62
std::vector< unsigned long > put() const
static double twoToMinus_53()
static bool checkFile(std::istream &file, const std::string &filename, const std::string &classname, const std::string &methodname)
Definition: RandomEngine.cc:45
static std::string beginTag()
ifstream in
Definition: comparison.C:7
static void getTheTableSeeds(long *seeds, int index)
Definition: Random.cc:251
void setSeeds(const long *seeds, int)
double flat()
Definition: G4AblaRandom.cc:48
void flatArray(const int size, double *vect)