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