CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: RandEngine.cc,v 1.8 2010/06/16 17:24:53 garren Exp $ 00002 // -*- C++ -*- 00003 // 00004 // ----------------------------------------------------------------------- 00005 // HEP Random 00006 // --- RandEngine --- 00007 // class implementation file 00008 // ----------------------------------------------------------------------- 00009 // This file is part of Geant4 (simulation toolkit for HEP). 00010 00011 // ======================================================================= 00012 // Gabriele Cosmo - Created: 5th September 1995 00013 // - Minor corrections: 31st October 1996 00014 // - Added methods for engine status: 19th November 1996 00015 // - mx is initialised to RAND_MAX: 2nd April 1997 00016 // - Fixed bug in setSeeds(): 15th September 1997 00017 // - Private copy constructor and operator=: 26th Feb 1998 00018 // J.Marraffino - Added stream operators and related constructor. 00019 // Added automatic seed selection from seed table and 00020 // engine counter. Removed RAND_MAX and replaced by 00021 // std::pow(0.5,32.). Flat() returns now 32 bits values 00022 // obtained by concatenation: 15th Feb 1998 00023 // Ken Smith - Added conversion operators: 6th Aug 1998 00024 // J. Marraffino - Remove dependence on hepString class 13 May 1999 00025 // M. Fischler - Rapaired bug that in flat() that relied on rand() to 00026 // deliver 15-bit results. This bug was causing flat() 00027 // on Linux systems to yield randoms with mean of 5/8(!) 00028 // - Modified algorithm such that on systems with 31-bit rand() 00029 // it will call rand() only once instead of twice. Feb 2004 00030 // M. Fischler - Modified the general-case template for RandEngineBuilder 00031 // such that when RAND_MAX is an unexpected value the routine 00032 // will still deliver a sensible flat() random. 00033 // M. Fischler - Methods for distrib. instance save/restore 12/8/04 00034 // M. Fischler - split get() into tag validation and 00035 // getState() for anonymous restores 12/27/04 00036 // M. Fischler - put/get for vectors of ulongs 3/14/05 00037 // M. Fischler - State-saving using only ints, for portability 4/12/05 00038 // 00039 // ======================================================================= 00040 00041 #include "CLHEP/Random/defs.h" 00042 #include "CLHEP/Random/RandEngine.h" 00043 #include "CLHEP/Random/Random.h" 00044 #include "CLHEP/Random/engineIDulong.h" 00045 #include <string.h> // for strcmp 00046 #include <cstdlib> // for int() 00047 00048 namespace CLHEP { 00049 00050 #ifdef NOTDEF 00051 // The way to test for proper behavior of the RandEngineBuilder 00052 // for arbitrary RAND_MAX, on a system where RAND_MAX is some 00053 // fixed specialized value and rand() behaves accordingly, is 00054 // to set up a fake RAND_MAX and a fake version of rand() 00055 // by enabling this block. 00056 #undef RAND_MAX 00057 #define RAND_MAX 9382956 00058 #include "CLHEP/Random/MTwistEngine.h" 00059 #include "CLHEP/Random/RandFlat.h" 00060 MTwistEngine * fakeFlat = new MTwistEngine; 00061 RandFlat rflat (fakeFlat, 0, RAND_MAX+1); 00062 int rand() { return (int)rflat(); } 00063 #endif 00064 00065 static const int MarkerLen = 64; // Enough room to hold a begin or end marker. 00066 00067 // number of instances with automatic seed selection 00068 int RandEngine::numEngines = 0; 00069 00070 // Maximum index into the seed table 00071 int RandEngine::maxIndex = 215; 00072 00073 std::string RandEngine::name() const {return "RandEngine";} 00074 00075 RandEngine::RandEngine(long seed) 00076 : HepRandomEngine() 00077 { 00078 setSeed(seed,0); 00079 setSeeds(&theSeed,0); 00080 seq = 0; 00081 } 00082 00083 RandEngine::RandEngine(int rowIndex, int colIndex) 00084 : HepRandomEngine() 00085 { 00086 long seeds[2]; 00087 long seed; 00088 00089 int cycle = std::abs(int(rowIndex/maxIndex)); 00090 int row = std::abs(int(rowIndex%maxIndex)); 00091 int col = std::abs(int(colIndex%2)); 00092 long mask = ((cycle & 0x000007ff) << 20 ); 00093 HepRandom::getTheTableSeeds( seeds, row ); 00094 seed = (seeds[col])^mask; 00095 setSeed(seed,0); 00096 setSeeds(&theSeed,0); 00097 seq = 0; 00098 } 00099 00100 RandEngine::RandEngine() 00101 : HepRandomEngine() 00102 { 00103 long seeds[2]; 00104 long seed; 00105 int cycle,curIndex; 00106 00107 cycle = std::abs(int(numEngines/maxIndex)); 00108 curIndex = std::abs(int(numEngines%maxIndex)); 00109 numEngines += 1; 00110 long mask = ((cycle & 0x007fffff) << 8); 00111 HepRandom::getTheTableSeeds( seeds, curIndex ); 00112 seed = seeds[0]^mask; 00113 setSeed(seed,0); 00114 setSeeds(&theSeed,0); 00115 seq = 0; 00116 } 00117 00118 RandEngine::RandEngine(std::istream& is) 00119 : HepRandomEngine() 00120 { 00121 is >> *this; 00122 } 00123 00124 RandEngine::~RandEngine() {} 00125 00126 void RandEngine::setSeed(long seed, int) 00127 { 00128 theSeed = seed; 00129 srand( int(seed) ); 00130 seq = 0; 00131 } 00132 00133 void RandEngine::setSeeds(const long* seeds, int) 00134 { 00135 setSeed(seeds ? *seeds : 19780503L, 0); 00136 theSeeds = seeds; 00137 } 00138 00139 void RandEngine::saveStatus( const char filename[] ) const 00140 { 00141 std::ofstream outFile( filename, std::ios::out ) ; 00142 00143 if (!outFile.bad()) { 00144 outFile << "Uvec\n"; 00145 std::vector<unsigned long> v = put(); 00146 #ifdef TRACE_IO 00147 std::cout << "Result of v = put() is:\n"; 00148 #endif 00149 for (unsigned int i=0; i<v.size(); ++i) { 00150 outFile << v[i] << "\n"; 00151 #ifdef TRACE_IO 00152 std::cout << v[i] << " "; 00153 if (i%6==0) std::cout << "\n"; 00154 #endif 00155 } 00156 #ifdef TRACE_IO 00157 std::cout << "\n"; 00158 #endif 00159 } 00160 #ifdef REMOVED 00161 if (!outFile.bad()) { 00162 outFile << theSeed << std::endl; 00163 outFile << seq << std::endl; 00164 } 00165 #endif 00166 } 00167 00168 void RandEngine::restoreStatus( const char filename[] ) 00169 { 00170 // The only way to restore the status of RandEngine is to 00171 // keep track of the number of shooted random sequences, reset 00172 // the engine and re-shoot them again. The Rand algorithm does 00173 // not provide any way of getting its internal status. 00174 00175 std::ifstream inFile( filename, std::ios::in); 00176 if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) { 00177 std::cout << " -- Engine state remains unchanged\n"; 00178 return; 00179 } 00180 if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) { 00181 std::vector<unsigned long> v; 00182 unsigned long xin; 00183 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) { 00184 inFile >> xin; 00185 #ifdef TRACE_IO 00186 std::cout << "ivec = " << ivec << " xin = " << xin << " "; 00187 if (ivec%3 == 0) std::cout << "\n"; 00188 #endif 00189 if (!inFile) { 00190 inFile.clear(std::ios::badbit | inFile.rdstate()); 00191 std::cerr << "\nRandEngine state (vector) description improper." 00192 << "\nrestoreStatus has failed." 00193 << "\nInput stream is probably mispositioned now." << std::endl; 00194 return; 00195 } 00196 v.push_back(xin); 00197 } 00198 getState(v); 00199 return; 00200 } 00201 00202 long count; 00203 00204 if (!inFile.bad() && !inFile.eof()) { 00205 // inFile >> theSeed; removed -- encompased by possibleKeywordInput 00206 inFile >> count; 00207 setSeed(theSeed,0); 00208 seq = 0; 00209 while (seq < count) flat(); 00210 } 00211 } 00212 00213 void RandEngine::showStatus() const 00214 { 00215 std::cout << std::endl; 00216 std::cout << "---------- Rand engine status ----------" << std::endl; 00217 std::cout << " Initial seed = " << theSeed << std::endl; 00218 std::cout << " Shooted sequences = " << seq << std::endl; 00219 std::cout << "----------------------------------------" << std::endl; 00220 } 00221 00222 // ==================================================== 00223 // Implementation of flat() (along with needed helpers) 00224 // ==================================================== 00225 00226 // Here we set up such that **at compile time**, the compiler decides based on 00227 // RAND_MAX how to form a random double with 32 leading random bits, using 00228 // one or two calls to rand(). Some of the lowest order bits of 32 are allowed 00229 // to be as weak as mere XORs of some higher bits, but not to be always fixed. 00230 // 00231 // The method decision is made at compile time, rather than using a run-time 00232 // if on the value of RAND_MAX. Although it is possible to cope with arbitrary 00233 // values of RAND_MAX of the form 2**N-1, with the same efficiency, the 00234 // template techniques needed would look mysterious and perhaps over-stress 00235 // some older compilers. We therefore only treat RAND_MAX = 2**15-1 (as on 00236 // most older systems) and 2**31-1 (as on the later Linux systems) as special 00237 // and super-efficient cases. We detect any different value, and use an 00238 // algorithm which is correct even if RAND_MAX is not one less than a power 00239 // of 2. 00240 00241 template <int> struct RandEngineBuilder { // RAND_MAX any arbitrary value 00242 static unsigned int thirtyTwoRandomBits(long& seq) { 00243 00244 static bool prepared = false; 00245 static unsigned int iT; 00246 static unsigned int iK; 00247 static unsigned int iS; 00248 static int iN; 00249 static double fS; 00250 static double fT; 00251 00252 if ( (RAND_MAX >> 31) > 0 ) 00253 { 00254 // Here, we know that integer arithmetic is 64 bits. 00255 if ( !prepared ) { 00256 iS = (unsigned long)RAND_MAX + 1; 00257 iK = 1; 00258 // int StoK = S; 00259 int StoK = iS; 00260 // The two statements below are equivalent, but some compilers 00261 // are getting too smart and complain about the first statement. 00262 //if ( (RAND_MAX >> 32) == 0) { 00263 if( (unsigned long) (RAND_MAX) <= (( (1uL) << 31 ) - 1 ) ) { 00264 iK = 2; 00265 // StoK = S*S; 00266 StoK = iS*iS; 00267 } 00268 int m; 00269 for ( m = 0; m < 64; ++m ) { 00270 StoK >>= 1; 00271 if (StoK == 0) break; 00272 } 00273 iT = 1 << m; 00274 prepared = true; 00275 } 00276 int v = 0; 00277 do { 00278 for ( int i = 0; i < iK; ++i ) { 00279 v = iS*v+rand(); ++seq; 00280 } 00281 } while (v < iT); 00282 return v & 0xFFFFFFFF; 00283 00284 } 00285 00286 else if ( (RAND_MAX >> 26) == 0 ) 00287 { 00288 // Here, we know it is safe to work in terms of doubles without loss 00289 // of precision, but we have no idea how many randoms we will need to 00290 // generate 32 bits. 00291 if ( !prepared ) { 00292 fS = (unsigned long)RAND_MAX + 1; 00293 double twoTo32 = std::ldexp(1.0,32); 00294 double StoK = fS; 00295 for ( iK = 1; StoK < twoTo32; StoK *= fS, iK++ ) { } 00296 int m; 00297 fT = 1.0; 00298 for ( m = 0; m < 64; ++m ) { 00299 StoK *= .5; 00300 if (StoK < 1.0) break; 00301 fT *= 2.0; 00302 } 00303 prepared = true; 00304 } 00305 double v = 0; 00306 do { 00307 for ( int i = 0; i < iK; ++i ) { 00308 v = fS*v+rand(); ++seq; 00309 } 00310 } while (v < fT); 00311 return ((unsigned int)v) & 0xFFFFFFFF; 00312 00313 } 00314 else 00315 { 00316 // Here, we know that 16 random bits are available from each of 00317 // two random numbers. 00318 if ( !prepared ) { 00319 iS = (unsigned long)RAND_MAX + 1; 00320 int SshiftN = iS; 00321 for (iN = 0; SshiftN > 1; SshiftN >>= 1, iN++) { } 00322 iN -= 17; 00323 prepared = true; 00324 } 00325 unsigned int x1, x2; 00326 do { x1 = rand(); ++seq;} while (x1 < (1<<16) ); 00327 do { x2 = rand(); ++seq;} while (x2 < (1<<16) ); 00328 return x1 | (x2 << 16); 00329 } 00330 00331 } 00332 }; 00333 00334 template <> struct RandEngineBuilder<2147483647> { // RAND_MAX = 2**31 - 1 00335 inline static unsigned int thirtyTwoRandomBits(long& seq) { 00336 unsigned int x = rand() << 1; ++seq; // bits 31-1 00337 x ^= ( (x>>23) ^ (x>>7) ) ^1; // bit 0 (weakly pseudo-random) 00338 return x & 0xFFFFFFFF; // mask in case int is 64 bits 00339 } 00340 }; 00341 00342 00343 template <> struct RandEngineBuilder<32767> { // RAND_MAX = 2**15 - 1 00344 inline static unsigned int thirtyTwoRandomBits(long& seq) { 00345 unsigned int x = rand() << 17; ++seq; // bits 31-17 00346 x ^= rand() << 2; ++seq; // bits 16-2 00347 x ^= ( (x>>23) ^ (x>>7) ) ^3; // bits 1-0 (weakly pseudo-random) 00348 return x & 0xFFFFFFFF; // mask in case int is 64 bits 00349 } 00350 }; 00351 00352 double RandEngine::flat() 00353 { 00354 double r; 00355 do { r = RandEngineBuilder<RAND_MAX>::thirtyTwoRandomBits(seq); 00356 } while ( r == 0 ); 00357 return r/4294967296.0; 00358 } 00359 00360 void RandEngine::flatArray(const int size, double* vect) 00361 { 00362 int i; 00363 00364 for (i=0; i<size; ++i) 00365 vect[i]=flat(); 00366 } 00367 00368 RandEngine::operator unsigned int() { 00369 return RandEngineBuilder<RAND_MAX>::thirtyTwoRandomBits(seq); 00370 } 00371 00372 std::ostream & RandEngine::put ( std::ostream& os ) const 00373 { 00374 char beginMarker[] = "RandEngine-begin"; 00375 char endMarker[] = "RandEngine-end"; 00376 00377 os << " " << beginMarker << "\n"; 00378 os << theSeed << " " << seq << " "; 00379 os << endMarker << "\n"; 00380 return os; 00381 } 00382 00383 std::vector<unsigned long> RandEngine::put () const { 00384 std::vector<unsigned long> v; 00385 v.push_back (engineIDulong<RandEngine>()); 00386 v.push_back(static_cast<unsigned long>(theSeed)); 00387 v.push_back(static_cast<unsigned long>(seq)); 00388 return v; 00389 } 00390 00391 std::istream & RandEngine::get ( std::istream& is ) 00392 { 00393 // The only way to restore the status of RandEngine is to 00394 // keep track of the number of shooted random sequences, reset 00395 // the engine and re-shoot them again. The Rand algorithm does 00396 // not provide any way of getting its internal status. 00397 char beginMarker [MarkerLen]; 00398 is >> std::ws; 00399 is.width(MarkerLen); // causes the next read to the char* to be <= 00400 // that many bytes, INCLUDING A TERMINATION \0 00401 // (Stroustrup, section 21.3.2) 00402 is >> beginMarker; 00403 if (strcmp(beginMarker,"RandEngine-begin")) { 00404 is.clear(std::ios::badbit | is.rdstate()); 00405 std::cout << "\nInput stream mispositioned or" 00406 << "\nRandEngine state description missing or" 00407 << "\nwrong engine type found." << std::endl; 00408 return is; 00409 } 00410 return getState(is); 00411 } 00412 00413 std::string RandEngine::beginTag ( ) { 00414 return "RandEngine-begin"; 00415 } 00416 00417 std::istream & RandEngine::getState ( std::istream& is ) 00418 { 00419 if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) { 00420 std::vector<unsigned long> v; 00421 unsigned long uu; 00422 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) { 00423 is >> uu; 00424 if (!is) { 00425 is.clear(std::ios::badbit | is.rdstate()); 00426 std::cerr << "\nRandEngine state (vector) description improper." 00427 << "\ngetState() has failed." 00428 << "\nInput stream is probably mispositioned now." << std::endl; 00429 return is; 00430 } 00431 v.push_back(uu); 00432 } 00433 getState(v); 00434 return (is); 00435 } 00436 00437 // is >> theSeed; Removed, encompassed by possibleKeywordInput() 00438 00439 char endMarker [MarkerLen]; 00440 long count; 00441 is >> count; 00442 is >> std::ws; 00443 is.width(MarkerLen); 00444 is >> endMarker; 00445 if (strcmp(endMarker,"RandEngine-end")) { 00446 is.clear(std::ios::badbit | is.rdstate()); 00447 std::cerr << "\nRandEngine state description incomplete." 00448 << "\nInput stream is probably mispositioned now." << std::endl; 00449 return is; 00450 } 00451 setSeed(theSeed,0); 00452 while (seq < count) flat(); 00453 return is; 00454 } 00455 00456 bool RandEngine::get (const std::vector<unsigned long> & v) { 00457 if ((v[0] & 0xffffffffUL) != engineIDulong<RandEngine>()) { 00458 std::cerr << 00459 "\nRandEngine get:state vector has wrong ID word - state unchanged\n"; 00460 return false; 00461 } 00462 return getState(v); 00463 } 00464 00465 bool RandEngine::getState (const std::vector<unsigned long> & v) { 00466 if (v.size() != VECTOR_STATE_SIZE ) { 00467 std::cerr << 00468 "\nRandEngine get:state vector has wrong length - state unchanged\n"; 00469 return false; 00470 } 00471 theSeed = v[1]; 00472 int count = v[2]; 00473 setSeed(theSeed,0); 00474 while (seq < count) flat(); 00475 return true; 00476 } 00477 } // namespace CLHEP