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