CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: MTwistEngine.cc,v 1.6 2010/06/16 17:24:53 garren Exp $ 00002 // -*- C++ -*- 00003 // 00004 // ----------------------------------------------------------------------- 00005 // HEP Random 00006 // --- MTwistEngine --- 00007 // class implementation file 00008 // ----------------------------------------------------------------------- 00009 // A "fast, compact, huge-period generator" based on M. Matsumoto and 00010 // T. Nishimura, "Mersenne Twister: A 623-dimensionally equidistributed 00011 // uniform pseudorandom number generator", to appear in ACM Trans. on 00012 // Modeling and Computer Simulation. It is a twisted GFSR generator 00013 // with a Mersenne-prime period of 2^19937-1, uniform on open interval (0,1) 00014 // ======================================================================= 00015 // Ken Smith - Started initial draft: 14th Jul 1998 00016 // - Optimized to get pow() out of flat() method 00017 // - Added conversion operators: 6th Aug 1998 00018 // J. Marraffino - Added some explicit casts to deal with 00019 // machines where sizeof(int) != sizeof(long) 22 Aug 1998 00020 // M. Fischler - Modified constructors such that no two 00021 // seeds will match sequences, no single/double 00022 // seeds will match, explicit seeds give 00023 // determined results, and default will not 00024 // match any of the above or other defaults. 00025 // - Modified use of the various exponents of 2 00026 // to avoid per-instance space overhead and 00027 // correct the rounding procedure 16 Sep 1998 00028 // J. Marfaffino - Remove dependence on hepString class 13 May 1999 00029 // M. Fischler - In restore, checkFile for file not found 03 Dec 2004 00030 // M. Fischler - Methods for distrib. instacne save/restore 12/8/04 00031 // M. Fischler - split get() into tag validation and 00032 // getState() for anonymous restores 12/27/04 00033 // M. Fischler - put/get for vectors of ulongs 3/14/05 00034 // M. Fischler - State-saving using only ints, for portability 4/12/05 00035 // M. Fischler - Improved seeding in setSeed (Savanah bug #17479) 11/15/06 00036 // - (Possible optimization - now that the seeding is improved, 00037 // is it still necessary for ctor to "warm up" by discarding 00038 // 2000 iterations?) 00039 // 00040 // ======================================================================= 00041 00042 #include "CLHEP/Random/defs.h" 00043 #include "CLHEP/Random/Random.h" 00044 #include "CLHEP/Random/MTwistEngine.h" 00045 #include "CLHEP/Random/engineIDulong.h" 00046 #include <string.h> // for strcmp 00047 #include <cstdlib> // for abs(int) 00048 00049 using namespace std; 00050 00051 namespace CLHEP { 00052 00053 static const int MarkerLen = 64; // Enough room to hold a begin or end marker. 00054 00055 std::string MTwistEngine::name() const {return "MTwistEngine";} 00056 00057 int MTwistEngine::numEngines = 0; 00058 int MTwistEngine::maxIndex = 215; 00059 00060 MTwistEngine::MTwistEngine() 00061 : HepRandomEngine() 00062 { 00063 int cycle = abs(int(numEngines/maxIndex)); 00064 int curIndex = abs(int(numEngines%maxIndex)); 00065 long mask = ((cycle & 0x007fffff) << 8); 00066 long seedlist[2]; 00067 HepRandom::getTheTableSeeds( seedlist, curIndex ); 00068 seedlist[0] = (seedlist[0])^mask; 00069 seedlist[1] = 0; 00070 setSeeds( seedlist, numEngines ); 00071 count624=0; 00072 ++numEngines; 00073 for( int i=0; i < 2000; ++i ) flat(); // Warm up just a bit 00074 } 00075 00076 MTwistEngine::MTwistEngine(long seed) 00077 : HepRandomEngine() 00078 { 00079 long seedlist[2]={seed,17587}; 00080 setSeeds( seedlist, 0 ); 00081 count624=0; 00082 for( int i=0; i < 2000; ++i ) flat(); // Warm up just a bit 00083 } 00084 00085 MTwistEngine::MTwistEngine(int rowIndex, int colIndex) 00086 : HepRandomEngine() 00087 { 00088 int cycle = abs(int(rowIndex/maxIndex)); 00089 int row = abs(int(rowIndex%maxIndex)); 00090 int col = abs(int(colIndex%2)); 00091 long mask = (( cycle & 0x000007ff ) << 20 ); 00092 long seedlist[2]; 00093 HepRandom::getTheTableSeeds( seedlist, row ); 00094 seedlist[0] = (seedlist[col])^mask; 00095 seedlist[1] = 690691; 00096 setSeeds(seedlist, 4444772); 00097 count624=0; 00098 for( int i=0; i < 2000; ++i ) flat(); // Warm up just a bit 00099 } 00100 00101 MTwistEngine::MTwistEngine( std::istream& is ) 00102 : HepRandomEngine() 00103 { 00104 is >> *this; 00105 } 00106 00107 MTwistEngine::~MTwistEngine() {} 00108 00109 double MTwistEngine::flat() { 00110 unsigned int y; 00111 00112 if( count624 >= N ) { 00113 register int i; 00114 00115 for( i=0; i < NminusM; ++i ) { 00116 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff); 00117 mt[i] = mt[i+M] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 ); 00118 } 00119 00120 for( ; i < N-1 ; ++i ) { 00121 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff); 00122 mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 ); 00123 } 00124 00125 y = (mt[i] & 0x80000000) | (mt[0] & 0x7fffffff); 00126 mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 ); 00127 00128 count624 = 0; 00129 } 00130 00131 y = mt[count624]; 00132 y ^= ( y >> 11); 00133 y ^= ((y << 7 ) & 0x9d2c5680); 00134 y ^= ((y << 15) & 0xefc60000); 00135 y ^= ( y >> 18); 00136 00137 return y * twoToMinus_32() + // Scale to range 00138 (mt[count624++] >> 11) * twoToMinus_53() + // fill remaining bits 00139 nearlyTwoToMinus_54(); // make sure non-zero 00140 } 00141 00142 void MTwistEngine::flatArray( const int size, double *vect ) { 00143 for( int i=0; i < size; ++i) vect[i] = flat(); 00144 } 00145 00146 void MTwistEngine::setSeed(long seed, int k) { 00147 00148 // MF 11/15/06 - Change seeding algorithm to a newer one recommended 00149 // by Matsumoto: The original algorithm was 00150 // mt[i] = (69069 * mt[i-1]) & 0xffffffff and this gives 00151 // problems when the seed bit pattern has lots of zeros 00152 // (for example, 0x08000000). Savanah bug #17479. 00153 00154 theSeed = seed ? seed : 4357; 00155 int mti; 00156 const int N1=624; 00157 mt[0] = (unsigned int) (theSeed&0xffffffffUL); 00158 for (mti=1; mti<N1; mti++) { 00159 mt[mti] = (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti); 00160 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ 00161 /* In the previous versions, MSBs of the seed affect */ 00162 /* only MSBs of the array mt[]. */ 00163 /* 2002/01/09 modified by Makoto Matsumoto */ 00164 mt[mti] &= 0xffffffffUL; 00165 /* for >32 bit machines */ 00166 } 00167 for( int i=1; i < 624; ++i ) { 00168 mt[i] ^= k; // MF 9/16/98: distinguish starting points 00169 } 00170 // MF 11/15/06 This distinction of starting points based on values of k 00171 // is kept even though the seeding algorithm itself is improved. 00172 } 00173 00174 void MTwistEngine::setSeeds(const long *seeds, int k) { 00175 setSeed( (*seeds ? *seeds : 43571346), k ); 00176 int i; 00177 for( i=1; i < 624; ++i ) { 00178 mt[i] = ( seeds[1] + mt[i] ) & 0xffffffff; // MF 9/16/98 00179 } 00180 theSeeds = seeds; 00181 } 00182 00183 void MTwistEngine::saveStatus( const char filename[] ) const 00184 { 00185 std::ofstream outFile( filename, std::ios::out ) ; 00186 if (!outFile.bad()) { 00187 outFile << theSeed << std::endl; 00188 for (int i=0; i<624; ++i) outFile <<std::setprecision(20) << mt[i] << " "; 00189 outFile << std::endl; 00190 outFile << count624 << std::endl; 00191 } 00192 } 00193 00194 void MTwistEngine::restoreStatus( const char filename[] ) 00195 { 00196 std::ifstream inFile( filename, std::ios::in); 00197 if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) { 00198 std::cerr << " -- Engine state remains unchanged\n"; 00199 return; 00200 } 00201 00202 if (!inFile.bad() && !inFile.eof()) { 00203 inFile >> theSeed; 00204 for (int i=0; i<624; ++i) inFile >> mt[i]; 00205 inFile >> count624; 00206 } 00207 } 00208 00209 void MTwistEngine::showStatus() const 00210 { 00211 std::cout << std::endl; 00212 std::cout << "--------- MTwist engine status ---------" << std::endl; 00213 std::cout << std::setprecision(20); 00214 std::cout << " Initial seed = " << theSeed << std::endl; 00215 std::cout << " Current index = " << count624 << std::endl; 00216 std::cout << " Array status mt[] = " << std::endl; 00217 for (int i=0; i<624; i+=5) { 00218 std::cout << mt[i] << " " << mt[i+1] << " " << mt[i+2] << " " 00219 << mt[i+3] << " " << mt[i+4] << std::endl; 00220 } 00221 std::cout << "----------------------------------------" << std::endl; 00222 } 00223 00224 MTwistEngine::operator float() { 00225 unsigned int y; 00226 00227 if( count624 >= N ) { 00228 register int i; 00229 00230 for( i=0; i < NminusM; ++i ) { 00231 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff); 00232 mt[i] = mt[i+M] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 ); 00233 } 00234 00235 for( ; i < N-1 ; ++i ) { 00236 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff); 00237 mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 ); 00238 } 00239 00240 y = (mt[i] & 0x80000000) | (mt[0] & 0x7fffffff); 00241 mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 ); 00242 00243 count624 = 0; 00244 } 00245 00246 y = mt[count624++]; 00247 y ^= ( y >> 11); 00248 y ^= ((y << 7 ) & 0x9d2c5680); 00249 y ^= ((y << 15) & 0xefc60000); 00250 y ^= ( y >> 18); 00251 00252 return (float)(y * twoToMinus_32()); 00253 } 00254 00255 MTwistEngine::operator unsigned int() { 00256 unsigned int y; 00257 00258 if( count624 >= N ) { 00259 register int i; 00260 00261 for( i=0; i < NminusM; ++i ) { 00262 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff); 00263 mt[i] = mt[i+M] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 ); 00264 } 00265 00266 for( ; i < N-1 ; ++i ) { 00267 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff); 00268 mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 ); 00269 } 00270 00271 y = (mt[i] & 0x80000000) | (mt[0] & 0x7fffffff); 00272 mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 ); 00273 00274 count624 = 0; 00275 } 00276 00277 y = mt[count624++]; 00278 y ^= ( y >> 11); 00279 y ^= ((y << 7 ) & 0x9d2c5680); 00280 y ^= ((y << 15) & 0xefc60000); 00281 y ^= ( y >> 18); 00282 00283 return y; 00284 } 00285 00286 std::ostream & MTwistEngine::put ( std::ostream& os ) const 00287 { 00288 char beginMarker[] = "MTwistEngine-begin"; 00289 char endMarker[] = "MTwistEngine-end"; 00290 00291 int pr = os.precision(20); 00292 os << " " << beginMarker << " "; 00293 os << theSeed << " "; 00294 for (int i=0; i<624; ++i) { 00295 os << mt[i] << "\n"; 00296 } 00297 os << count624 << " "; 00298 os << endMarker << "\n"; 00299 os.precision(pr); 00300 return os; 00301 } 00302 00303 std::vector<unsigned long> MTwistEngine::put () const { 00304 std::vector<unsigned long> v; 00305 v.push_back (engineIDulong<MTwistEngine>()); 00306 for (int i=0; i<624; ++i) { 00307 v.push_back(static_cast<unsigned long>(mt[i])); 00308 } 00309 v.push_back(count624); 00310 return v; 00311 } 00312 00313 std::istream & MTwistEngine::get ( std::istream& is ) 00314 { 00315 char beginMarker [MarkerLen]; 00316 is >> std::ws; 00317 is.width(MarkerLen); // causes the next read to the char* to be <= 00318 // that many bytes, INCLUDING A TERMINATION \0 00319 // (Stroustrup, section 21.3.2) 00320 is >> beginMarker; 00321 if (strcmp(beginMarker,"MTwistEngine-begin")) { 00322 is.clear(std::ios::badbit | is.rdstate()); 00323 std::cerr << "\nInput stream mispositioned or" 00324 << "\nMTwistEngine state description missing or" 00325 << "\nwrong engine type found." << std::endl; 00326 return is; 00327 } 00328 return getState(is); 00329 } 00330 00331 std::string MTwistEngine::beginTag ( ) { 00332 return "MTwistEngine-begin"; 00333 } 00334 00335 std::istream & MTwistEngine::getState ( std::istream& is ) 00336 { 00337 char endMarker [MarkerLen]; 00338 is >> theSeed; 00339 for (int i=0; i<624; ++i) is >> mt[i]; 00340 is >> count624; 00341 is >> std::ws; 00342 is.width(MarkerLen); 00343 is >> endMarker; 00344 if (strcmp(endMarker,"MTwistEngine-end")) { 00345 is.clear(std::ios::badbit | is.rdstate()); 00346 std::cerr << "\nMTwistEngine state description incomplete." 00347 << "\nInput stream is probably mispositioned now." << std::endl; 00348 return is; 00349 } 00350 return is; 00351 } 00352 00353 bool MTwistEngine::get (const std::vector<unsigned long> & v) { 00354 if ((v[0] & 0xffffffffUL) != engineIDulong<MTwistEngine>()) { 00355 std::cerr << 00356 "\nMTwistEngine get:state vector has wrong ID word - state unchanged\n"; 00357 return false; 00358 } 00359 return getState(v); 00360 } 00361 00362 bool MTwistEngine::getState (const std::vector<unsigned long> & v) { 00363 if (v.size() != VECTOR_STATE_SIZE ) { 00364 std::cerr << 00365 "\nMTwistEngine get:state vector has wrong length - state unchanged\n"; 00366 return false; 00367 } 00368 for (int i=0; i<624; ++i) { 00369 mt[i]=v[i+1]; 00370 } 00371 count624 = v[625]; 00372 return true; 00373 } 00374 00375 } // namespace CLHEP