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