| Geant4 Cross Reference |
1 // -*- C++ -*- 1
2 //
3 // -------------------------------------------
4 // HEP Random
5 // --- MTwistEngine ---
6 // class implementation f
7 // -------------------------------------------
8 // A "fast, compact, huge-period generator" ba
9 // T. Nishimura, "Mersenne Twister: A 623-dime
10 // uniform pseudorandom number generator", to
11 // Modeling and Computer Simulation. It is a
12 // with a Mersenne-prime period of 2^19937-1,
13 // ===========================================
14 // Ken Smith - Started initial draft:
15 // - Optimized to get std::pow(
16 // - Added conversion operators
17 // J. Marraffino - Added some explicit casts
18 // machines where sizeof(int)
19 // M. Fischler - Modified constructors such
20 // seeds will match sequences, no singl
21 // seeds will match, explicit seeds giv
22 // determined results, and default will
23 // match any of the above or other defa
24 // - Modified use of the variou
25 // to avoid per-instance spac
26 // correct the rounding proce
27 // J. Marraffino - Remove dependence on hepSt
28 // M. Fischler - In restore, checkFile for
29 // M. Fischler - Methods for distrib. insta
30 // M. Fischler - split get() into tag valid
31 // getState() for anonymous r
32 // M. Fischler - put/get for vectors of ulo
33 // M. Fischler - State-saving using only in
34 // M. Fischler - Improved seeding in setSee
35 // - (Possible optimization - now that th
36 // is it still necessary for ctor to "w
37 // 2000 iterations?)
38 //
39 // ===========================================
40
41 #include "CLHEP/Random/Random.h"
42 #include "CLHEP/Random/MTwistEngine.h"
43 #include "CLHEP/Random/engineIDulong.h"
44 #include "CLHEP/Utility/atomic_int.h"
45
46 #include <atomic>
47 #include <cmath>
48 #include <iostream>
49 #include <string.h> // for strcmp
50 #include <vector>
51
52 namespace CLHEP {
53
54 namespace {
55 // Number of instances with automatic seed s
56 CLHEP_ATOMIC_INT_TYPE numberOfEngines(0);
57
58 // Maximum index into the seed table
59 const int maxIndex = 215;
60 }
61
62 static const int MarkerLen = 64; // Enough roo
63
64 std::string MTwistEngine::name() const {return
65
66 MTwistEngine::MTwistEngine()
67 : HepRandomEngine()
68 {
69 int numEngines = numberOfEngines++;
70 int cycle = std::abs(int(numEngines/maxIndex
71 int curIndex = std::abs(int(numEngines%maxIn
72 long mask = ((cycle & 0x007fffff) << 8);
73 long seedlist[2];
74 HepRandom::getTheTableSeeds( seedlist, curIn
75 seedlist[0] = (seedlist[0])^mask;
76 seedlist[1] = 0;
77 setSeeds( seedlist, numEngines );
78 count624=0;
79
80 for( int i=0; i < 2000; ++i ) flat(); /
81 }
82
83 MTwistEngine::MTwistEngine(long seed)
84 : HepRandomEngine()
85 {
86 long seedlist[2]={seed,17587};
87 setSeeds( seedlist, 0 );
88 count624=0;
89 for( int i=0; i < 2000; ++i ) flat(); /
90 }
91
92 MTwistEngine::MTwistEngine(int rowIndex, int c
93 : HepRandomEngine()
94 {
95 int cycle = std::abs(int(rowIndex/maxIndex))
96 int row = std::abs(int(rowIndex%maxIndex));
97 int col = std::abs(int(colIndex%2));
98 long mask = (( cycle & 0x000007ff ) << 20 );
99 long seedlist[2];
100 HepRandom::getTheTableSeeds( seedlist, row )
101 seedlist[0] = (seedlist[col])^mask;
102 seedlist[1] = 690691;
103 setSeeds(seedlist, 4444772);
104 count624=0;
105 for( int i=0; i < 2000; ++i ) flat(); /
106 }
107
108 MTwistEngine::MTwistEngine( std::istream& is )
109 : HepRandomEngine()
110 {
111 is >> *this;
112 }
113
114 MTwistEngine::~MTwistEngine() {}
115
116 double MTwistEngine::flat() {
117 unsigned int y;
118
119 if( count624 >= N ) {
120 int i;
121
122 for( i=0; i < NminusM; ++i ) {
123 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x
124 mt[i] = mt[i+M] ^ (y >> 1) ^ ((y &
125 }
126
127 for( ; i < N-1 ; ++i ) {
128 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x
129 mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y &
130 }
131
132 y = (mt[i] & 0x80000000) | (mt[0] & 0x7fff
133 mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ?
134
135 count624 = 0;
136 }
137
138 y = mt[count624];
139 y ^= ( y >> 11);
140 y ^= ((y << 7 ) & 0x9d2c5680);
141 y ^= ((y << 15) & 0xefc60000);
142 y ^= ( y >> 18);
143
144 return y * twoToMinus_3
145 (mt[count624++] >> 11) * twoToMinus_5
146 nearlyTwoToMinus_54();
147 }
148
149 void MTwistEngine::flatArray( const int size,
150 for( int i=0; i < size; ++i) vect[i] = flat(
151 }
152
153 void MTwistEngine::setSeed(long seed, int k) {
154
155 // MF 11/15/06 - Change seeding algorithm to
156 // by Matsumoto: The original
157 // mt[i] = (69069 * mt[i-1]) & 0xfffff
158 // problems when the seed bit pattern
159 // (for example, 0x08000000). Savanah
160
161 theSeed = seed ? seed : 4357;
162 int mti;
163 const int N1=624;
164 mt[0] = (unsigned int) (theSeed&0xffffffffUL
165 for (mti=1; mti<N1; mti++) {
166 mt[mti] = (1812433253UL * (mt[mti-1] ^ (mt
167 /* See Knuth TAOCP Vol2. 3rd Ed. P.106
168 /* In the previous versions, MSBs of t
169 /* only MSBs of the array mt[].
170 /* 2002/01/09 modified by Makoto Matsu
171 mt[mti] &= 0xffffffffUL;
172 /* for >32 bit machines */
173 }
174 for( int i=1; i < 624; ++i ) {
175 mt[i] ^= k; // MF 9/16/98: distinguish
176 }
177 // MF 11/15/06 This distinction of starting
178 // is kept even though the seedi
179 }
180
181 void MTwistEngine::setSeeds(const long *seeds,
182 setSeed( (*seeds ? *seeds : 43571346), k );
183 int i;
184 for( i=1; i < 624; ++i ) {
185 mt[i] = ( seeds[1] + mt[i] ) & 0xffffffff;
186 }
187 theSeeds = seeds;
188 }
189
190 void MTwistEngine::saveStatus( const char file
191 {
192 std::ofstream outFile( filename, std::ios::
193 if (!outFile.bad()) {
194 outFile << theSeed << std::endl;
195 for (int i=0; i<624; ++i) outFile <<std::
196 outFile << std::endl;
197 outFile << count624 << std::endl;
198 }
199 }
200
201 void MTwistEngine::restoreStatus( const char f
202 {
203 std::ifstream inFile( filename, std::ios::i
204 if (!checkFile ( inFile, filename, engineNa
205 std::cerr << " -- Engine state remains u
206 return;
207 }
208
209 if (!inFile.bad() && !inFile.eof()) {
210 inFile >> theSeed;
211 for (int i=0; i<624; ++i) inFile >> mt[i]
212 inFile >> count624;
213 }
214 }
215
216 void MTwistEngine::showStatus() const
217 {
218 std::cout << std::endl;
219 std::cout << "--------- MTwist engine statu
220 std::cout << std::setprecision(20);
221 std::cout << " Initial seed = " << the
222 std::cout << " Current index = " << cou
223 std::cout << " Array status mt[] = " << std
224 // 2014/06/06 L Garren
225 // the final line has 4 elements, not 5
226 for (int i=0; i<620; i+=5) {
227 std::cout << mt[i] << " " << mt[i+1] <<
228 << mt[i+3] << " " << mt[i+4] << "\n";
229 }
230 std::cout << mt[620] << " " << mt[621] <<
231 << mt[623] << std::endl;
232 std::cout << "-----------------------------
233 }
234
235 MTwistEngine::operator double() {
236 return flat();
237 }
238
239 MTwistEngine::operator float() {
240 unsigned int y;
241
242 if( count624 >= N ) {
243 int i;
244
245 for( i=0; i < NminusM; ++i ) {
246 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x
247 mt[i] = mt[i+M] ^ (y >> 1) ^ ((y &
248 }
249
250 for( ; i < N-1 ; ++i ) {
251 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x
252 mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y &
253 }
254
255 y = (mt[i] & 0x80000000) | (mt[0] & 0x7fff
256 mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ?
257
258 count624 = 0;
259 }
260
261 y = mt[count624++];
262 y ^= ( y >> 11);
263 y ^= ((y << 7 ) & 0x9d2c5680);
264 y ^= ((y << 15) & 0xefc60000);
265 y ^= ( y >> 18);
266
267 return (float)(y * twoToMinus_32());
268 }
269
270 MTwistEngine::operator unsigned int() {
271 unsigned int y;
272
273 if( count624 >= N ) {
274 int i;
275
276 for( i=0; i < NminusM; ++i ) {
277 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x
278 mt[i] = mt[i+M] ^ (y >> 1) ^ ((y &
279 }
280
281 for( ; i < N-1 ; ++i ) {
282 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x
283 mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y &
284 }
285
286 y = (mt[i] & 0x80000000) | (mt[0] & 0x7fff
287 mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ?
288
289 count624 = 0;
290 }
291
292 y = mt[count624++];
293 y ^= ( y >> 11);
294 y ^= ((y << 7 ) & 0x9d2c5680);
295 y ^= ((y << 15) & 0xefc60000);
296 y ^= ( y >> 18);
297
298 return y;
299 }
300
301 std::ostream & MTwistEngine::put ( std::ostrea
302 {
303 char beginMarker[] = "MTwistEngine-begin";
304 char endMarker[] = "MTwistEngine-end";
305
306 long pr = os.precision(20);
307 os << " " << beginMarker << " ";
308 os << theSeed << " ";
309 for (int i=0; i<624; ++i) {
310 os << mt[i] << "\n";
311 }
312 os << count624 << " ";
313 os << endMarker << "\n";
314 os.precision(pr);
315 return os;
316 }
317
318 std::vector<unsigned long> MTwistEngine::put (
319 std::vector<unsigned long> v;
320 v.push_back (engineIDulong<MTwistEngine>());
321 for (int i=0; i<624; ++i) {
322 v.push_back(static_cast<unsigned long>(mt
323 }
324 v.push_back(count624);
325 return v;
326 }
327
328 std::istream & MTwistEngine::get ( std::istre
329 {
330 char beginMarker [MarkerLen];
331 is >> std::ws;
332 is.width(MarkerLen); // causes the next rea
333 // that many bytes, INCLUDING A TERMINAT
334 // (Stroustrup, section 21.3.2)
335 is >> beginMarker;
336 if (strcmp(beginMarker,"MTwistEngine-begin")
337 is.clear(std::ios::badbit | is.rdstate())
338 std::cerr << "\nInput stream mispositione
339 << "\nMTwistEngine state description
340 << "\nwrong engine type found." << st
341 return is;
342 }
343 return getState(is);
344 }
345
346 std::string MTwistEngine::beginTag ( ) {
347 return "MTwistEngine-begin";
348 }
349
350 std::istream & MTwistEngine::getState ( std::
351 {
352 char endMarker [MarkerLen];
353 is >> theSeed;
354 for (int i=0; i<624; ++i) is >> mt[i];
355 is >> count624;
356 is >> std::ws;
357 is.width(MarkerLen);
358 is >> endMarker;
359 if (strcmp(endMarker,"MTwistEngine-end")) {
360 is.clear(std::ios::badbit | is.rdstate())
361 std::cerr << "\nMTwistEngine state descri
362 << "\nInput stream is probably mispos
363 return is;
364 }
365 return is;
366 }
367
368 bool MTwistEngine::get (const std::vector<unsi
369 if ((v[0] & 0xffffffffUL) != engineIDulong<M
370 std::cerr <<
371 "\nMTwistEngine get:state vector has wro
372 return false;
373 }
374 return getState(v);
375 }
376
377 bool MTwistEngine::getState (const std::vector
378 if (v.size() != VECTOR_STATE_SIZE ) {
379 std::cerr <<
380 "\nMTwistEngine get:state vector has wro
381 return false;
382 }
383 for (int i=0; i<624; ++i) {
384 mt[i]=(unsigned int)v[i+1];
385 }
386 count624 = (int)v[625];
387 return true;
388 }
389
390 } // namespace CLHEP
391