CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: MTwistEngine.cc,v 1.4.4.3.2.1 2008/11/13 21:35:23 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> 00047 #include <cmath> // for ldexp() 00048 #include <stdlib.h> // for abs(int) 00049 00050 using namespace std; 00051 00052 namespace CLHEP { 00053 00054 static const int MarkerLen = 64; // Enough room to hold a begin or end marker. 00055 00056 double MTwistEngine::twoToMinus_32; 00057 double MTwistEngine::twoToMinus_53; 00058 double MTwistEngine::nearlyTwoToMinus_54; 00059 00060 std::string MTwistEngine::name() const {return "MTwistEngine";} 00061 00062 00063 void MTwistEngine::powersOfTwo() { 00064 twoToMinus_32 = ldexp (1.0, -32); 00065 twoToMinus_53 = ldexp (1.0, -53); 00066 nearlyTwoToMinus_54 = ldexp (1.0, -54) - ldexp (1.0, -100); 00067 } 00068 00069 int MTwistEngine::numEngines = 0; 00070 int MTwistEngine::maxIndex = 215; 00071 00072 MTwistEngine::MTwistEngine() 00073 { 00074 powersOfTwo(); 00075 int cycle = abs(int(numEngines/maxIndex)); 00076 int curIndex = abs(int(numEngines%maxIndex)); 00077 long mask = ((cycle & 0x007fffff) << 8); 00078 long seedlist[2]; 00079 HepRandom::getTheTableSeeds( seedlist, curIndex ); 00080 seedlist[0] = (seedlist[0])^mask; 00081 seedlist[1] = 0; 00082 setSeeds( seedlist, numEngines ); 00083 count624=0; 00084 ++numEngines; 00085 for( int i=0; i < 2000; ++i ) flat(); // Warm up just a bit 00086 } 00087 00088 MTwistEngine::MTwistEngine(long seed) 00089 { 00090 powersOfTwo(); 00091 long seedlist[2]={seed,17587}; 00092 setSeeds( seedlist, 0 ); 00093 count624=0; 00094 for( int i=0; i < 2000; ++i ) flat(); // Warm up just a bit 00095 } 00096 00097 MTwistEngine::MTwistEngine(int rowIndex, int colIndex) 00098 { 00099 powersOfTwo(); 00100 int cycle = abs(int(rowIndex/maxIndex)); 00101 int row = abs(int(rowIndex%maxIndex)); 00102 int col = abs(int(colIndex%2)); 00103 long mask = (( cycle & 0x000007ff ) << 20 ); 00104 long seedlist[2]; 00105 HepRandom::getTheTableSeeds( seedlist, row ); 00106 seedlist[0] = (seedlist[col])^mask; 00107 seedlist[1] = 690691; 00108 setSeeds(seedlist, 4444772); 00109 count624=0; 00110 for( int i=0; i < 2000; ++i ) flat(); // Warm up just a bit 00111 } 00112 00113 MTwistEngine::MTwistEngine( std::istream& is ) 00114 { 00115 is >> *this; 00116 } 00117 00118 MTwistEngine::~MTwistEngine() {} 00119 00120 MTwistEngine::MTwistEngine( const MTwistEngine & p ) 00121 { 00122 *this = p; 00123 } 00124 00125 MTwistEngine & MTwistEngine::operator=( const MTwistEngine & p ) { 00126 if( this != &p ) { 00127 for( int i=0; i < 624; ++i ) mt[i] = p.mt[i]; 00128 count624 = p.count624; 00129 } 00130 return *this; 00131 } 00132 00133 double MTwistEngine::flat() { 00134 unsigned int y; 00135 00136 if( count624 >= N ) { 00137 register int i; 00138 00139 for( i=0; i < NminusM; ++i ) { 00140 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff); 00141 mt[i] = mt[i+M] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 ); 00142 } 00143 00144 for( ; i < N-1 ; ++i ) { 00145 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff); 00146 mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 ); 00147 } 00148 00149 y = (mt[i] & 0x80000000) | (mt[0] & 0x7fffffff); 00150 mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 ); 00151 00152 count624 = 0; 00153 } 00154 00155 y = mt[count624]; 00156 y ^= ( y >> 11); 00157 y ^= ((y << 7 ) & 0x9d2c5680); 00158 y ^= ((y << 15) & 0xefc60000); 00159 y ^= ( y >> 18); 00160 00161 return y * twoToMinus_32 + // Scale to range 00162 (mt[count624++] >> 11) * twoToMinus_53 + // fill remaining bits 00163 nearlyTwoToMinus_54; // make sure non-zero 00164 } 00165 00166 void MTwistEngine::flatArray( const int size, double *vect ) { 00167 for( int i=0; i < size; ++i) vect[i] = flat(); 00168 } 00169 00170 void MTwistEngine::setSeed(long seed, int k) { 00171 00172 // MF 11/15/06 - Change seeding algorithm to a newer one recommended 00173 // by Matsumoto: The original algorithm was 00174 // mt[i] = (69069 * mt[i-1]) & 0xffffffff and this gives 00175 // problems when the seed bit pattern has lots of zeros 00176 // (for example, 0x08000000). Savanah bug #17479. 00177 00178 theSeed = seed ? seed : 4357; 00179 int mti; 00180 const int N=624; 00181 mt[0] = (unsigned int) (theSeed&0xffffffffUL); 00182 for (mti=1; mti<N; mti++) { 00183 mt[mti] = (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti); 00184 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ 00185 /* In the previous versions, MSBs of the seed affect */ 00186 /* only MSBs of the array mt[]. */ 00187 /* 2002/01/09 modified by Makoto Matsumoto */ 00188 mt[mti] &= 0xffffffffUL; 00189 /* for >32 bit machines */ 00190 } 00191 for( int i=1; i < 624; ++i ) { 00192 mt[i] ^= k; // MF 9/16/98: distinguish starting points 00193 } 00194 // MF 11/15/06 This distinction of starting points based on values of k 00195 // is kept even though the seeding algorithm itself is improved. 00196 } 00197 00198 void MTwistEngine::setSeeds(const long *seeds, int k) { 00199 setSeed( (*seeds ? *seeds : 43571346), k ); 00200 int i; 00201 for( i=1; i < 624; ++i ) { 00202 mt[i] = ( seeds[1] + mt[i] ) & 0xffffffff; // MF 9/16/98 00203 } 00204 theSeeds = seeds; 00205 } 00206 00207 void MTwistEngine::saveStatus( const char filename[] ) const 00208 { 00209 std::ofstream outFile( filename, std::ios::out ) ; 00210 if (!outFile.bad()) { 00211 outFile << theSeed << std::endl; 00212 for (int i=0; i<624; ++i) outFile <<std::setprecision(20) << mt[i] << " "; 00213 outFile << std::endl; 00214 outFile << count624 << std::endl; 00215 } 00216 } 00217 00218 void MTwistEngine::restoreStatus( const char filename[] ) 00219 { 00220 std::ifstream inFile( filename, std::ios::in); 00221 if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) { 00222 std::cerr << " -- Engine state remains unchanged\n"; 00223 return; 00224 } 00225 00226 if (!inFile.bad() && !inFile.eof()) { 00227 inFile >> theSeed; 00228 for (int i=0; i<624; ++i) inFile >> mt[i]; 00229 inFile >> count624; 00230 } 00231 } 00232 00233 void MTwistEngine::showStatus() const 00234 { 00235 std::cout << std::endl; 00236 std::cout << "--------- MTwist engine status ---------" << std::endl; 00237 std::cout << std::setprecision(20); 00238 std::cout << " Initial seed = " << theSeed << std::endl; 00239 std::cout << " Current index = " << count624 << std::endl; 00240 std::cout << " Array status mt[] = " << std::endl; 00241 for (int i=0; i<624; i+=5) { 00242 std::cout << mt[i] << " " << mt[i+1] << " " << mt[i+2] << " " 00243 << mt[i+3] << " " << mt[i+4] << std::endl; 00244 } 00245 std::cout << "----------------------------------------" << std::endl; 00246 } 00247 00248 MTwistEngine::operator float() { 00249 unsigned int y; 00250 00251 if( count624 >= N ) { 00252 register int i; 00253 00254 for( i=0; i < NminusM; ++i ) { 00255 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff); 00256 mt[i] = mt[i+M] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 ); 00257 } 00258 00259 for( ; i < N-1 ; ++i ) { 00260 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff); 00261 mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 ); 00262 } 00263 00264 y = (mt[i] & 0x80000000) | (mt[0] & 0x7fffffff); 00265 mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 ); 00266 00267 count624 = 0; 00268 } 00269 00270 y = mt[count624++]; 00271 y ^= ( y >> 11); 00272 y ^= ((y << 7 ) & 0x9d2c5680); 00273 y ^= ((y << 15) & 0xefc60000); 00274 y ^= ( y >> 18); 00275 00276 return (float)(y * twoToMinus_32); 00277 } 00278 00279 MTwistEngine::operator unsigned int() { 00280 unsigned int y; 00281 00282 if( count624 >= N ) { 00283 register int i; 00284 00285 for( i=0; i < NminusM; ++i ) { 00286 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff); 00287 mt[i] = mt[i+M] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 ); 00288 } 00289 00290 for( ; i < N-1 ; ++i ) { 00291 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff); 00292 mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 ); 00293 } 00294 00295 y = (mt[i] & 0x80000000) | (mt[0] & 0x7fffffff); 00296 mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 ); 00297 00298 count624 = 0; 00299 } 00300 00301 y = mt[count624++]; 00302 y ^= ( y >> 11); 00303 y ^= ((y << 7 ) & 0x9d2c5680); 00304 y ^= ((y << 15) & 0xefc60000); 00305 y ^= ( y >> 18); 00306 00307 return y; 00308 } 00309 00310 std::ostream & MTwistEngine::put ( std::ostream& os ) const 00311 { 00312 char beginMarker[] = "MTwistEngine-begin"; 00313 char endMarker[] = "MTwistEngine-end"; 00314 00315 int pr = os.precision(20); 00316 os << " " << beginMarker << " "; 00317 os << theSeed << " "; 00318 for (int i=0; i<624; ++i) { 00319 os << mt[i] << "\n"; 00320 } 00321 os << count624 << " "; 00322 os << endMarker << "\n"; 00323 os.precision(pr); 00324 return os; 00325 } 00326 00327 std::vector<unsigned long> MTwistEngine::put () const { 00328 std::vector<unsigned long> v; 00329 v.push_back (engineIDulong<MTwistEngine>()); 00330 for (int i=0; i<624; ++i) { 00331 v.push_back(static_cast<unsigned long>(mt[i])); 00332 } 00333 v.push_back(count624); 00334 return v; 00335 } 00336 00337 std::istream & MTwistEngine::get ( std::istream& is ) 00338 { 00339 char beginMarker [MarkerLen]; 00340 is >> std::ws; 00341 is.width(MarkerLen); // causes the next read to the char* to be <= 00342 // that many bytes, INCLUDING A TERMINATION \0 00343 // (Stroustrup, section 21.3.2) 00344 is >> beginMarker; 00345 if (strcmp(beginMarker,"MTwistEngine-begin")) { 00346 is.clear(std::ios::badbit | is.rdstate()); 00347 std::cerr << "\nInput stream mispositioned or" 00348 << "\nMTwistEngine state description missing or" 00349 << "\nwrong engine type found." << std::endl; 00350 return is; 00351 } 00352 return getState(is); 00353 } 00354 00355 std::string MTwistEngine::beginTag ( ) { 00356 return "MTwistEngine-begin"; 00357 } 00358 00359 std::istream & MTwistEngine::getState ( std::istream& is ) 00360 { 00361 char endMarker [MarkerLen]; 00362 is >> theSeed; 00363 for (int i=0; i<624; ++i) is >> mt[i]; 00364 is >> count624; 00365 is >> std::ws; 00366 is.width(MarkerLen); 00367 is >> endMarker; 00368 if (strcmp(endMarker,"MTwistEngine-end")) { 00369 is.clear(std::ios::badbit | is.rdstate()); 00370 std::cerr << "\nMTwistEngine state description incomplete." 00371 << "\nInput stream is probably mispositioned now." << std::endl; 00372 return is; 00373 } 00374 return is; 00375 } 00376 00377 bool MTwistEngine::get (const std::vector<unsigned long> & v) { 00378 if ((v[0] & 0xffffffffUL) != engineIDulong<MTwistEngine>()) { 00379 std::cerr << 00380 "\nMTwistEngine get:state vector has wrong ID word - state unchanged\n"; 00381 return false; 00382 } 00383 return getState(v); 00384 } 00385 00386 bool MTwistEngine::getState (const std::vector<unsigned long> & v) { 00387 if (v.size() != VECTOR_STATE_SIZE ) { 00388 std::cerr << 00389 "\nMTwistEngine get:state vector has wrong length - state unchanged\n"; 00390 return false; 00391 } 00392 for (int i=0; i<624; ++i) { 00393 mt[i]=v[i+1]; 00394 } 00395 count624 = v[625]; 00396 return true; 00397 } 00398 00399 } // namespace CLHEP