CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

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