CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: TripleRand.cc,v 1.6 2010/06/16 17:24:53 garren Exp $ 00002 // -*- C++ -*- 00003 // 00004 // ----------------------------------------------------------------------- 00005 // Hep Random 00006 // --- TripleRand --- 00007 // class implementation file 00008 // ----------------------------------------------------------------------- 00009 // A canopy pseudo-random number generator. Using the Tausworthe 00010 // exclusive-or shift register, a simple Integer Coungruence generator, and 00011 // the Hurd 288 total bit shift register, all XOR'd with each other, we 00012 // provide an engine that should be a fairly good "mother" generator. 00013 // ======================================================================= 00014 // Ken Smith - Initial draft started: 23rd Jul 1998 00015 // - Added conversion operators: 6th Aug 1998 00016 // J. Marraffino - Added some explicit casts to deal with 00017 // machines where sizeof(int) != sizeof(long) 22 Aug 1998 00018 // M. Fischler - Modified use of the various exponents of 2 00019 // to avoid per-instance space overhead and 00020 // correct the rounding procedure 15 Sep 1998 00021 // - modify constructors so that no sequence can 00022 // ever accidentally be produced by differnt 00023 // seeds, and so that exxcept for the default 00024 // constructor, the seed fully determines the 00025 // sequence. 15 Sep 1998 00026 // M. Fischler - Eliminated dependency on hepString 13 May 1999 00027 // M. Fischler - Eliminated Taus() and Cong() which were 00028 // causing spurions warnings on SUN CC 27 May 1999 00029 // M. Fischler - Put endl at end of puts of Tausworth and 10 Apr 2001 00030 // integerCong 00031 // M. Fischler - In restore, checkFile for file not found 03 Dec 2004 00032 // M. Fischler - Methods put, get for instance save/restore 12/8/04 00033 // M. Fischler - split get() into tag validation and 00034 // getState() for anonymous restores 12/27/04 00035 // M. Fischler - put/get for vectors of ulongs 3/14/05 00036 // M. Fischler - State-saving using only ints, for portability 4/12/05 00037 // 00038 // ======================================================================= 00039 00040 #include "CLHEP/Random/TripleRand.h" 00041 #include "CLHEP/Random/defs.h" 00042 #include "CLHEP/Random/engineIDulong.h" 00043 #include <string.h> // for strcmp 00044 00045 namespace CLHEP { 00046 00047 static const int MarkerLen = 64; // Enough room to hold a begin or end marker. 00048 00049 //******************************************************************** 00050 // TripleRand 00051 //******************************************************************** 00052 00053 // Number of instances with automatic seed selection 00054 int TripleRand::numEngines = 0; 00055 00056 std::string TripleRand::name() const {return "TripleRand";} 00057 00058 TripleRand::TripleRand() 00059 : HepRandomEngine(), 00060 tausworthe (1234567 + numEngines + 175321), 00061 integerCong(69607 * tausworthe + 54329, numEngines), 00062 hurd(19781127 + integerCong) 00063 { 00064 theSeed = 1234567; 00065 ++numEngines; 00066 } 00067 00068 TripleRand::TripleRand(long seed) 00069 : HepRandomEngine(), 00070 tausworthe ((unsigned int)seed + 175321), 00071 integerCong(69607 * tausworthe + 54329, 1313), 00072 hurd(19781127 + integerCong) 00073 { 00074 theSeed = seed; 00075 } 00076 00077 TripleRand::TripleRand(std::istream & is) 00078 : HepRandomEngine() 00079 { 00080 is >> *this; 00081 } 00082 00083 TripleRand::TripleRand(int rowIndex, int colIndex) 00084 : HepRandomEngine(), 00085 tausworthe (rowIndex + numEngines * colIndex + 175321), 00086 integerCong(69607 * tausworthe + 54329, 19), 00087 hurd(19781127 + integerCong) 00088 { 00089 theSeed = rowIndex; 00090 } 00091 00092 TripleRand::~TripleRand() { } 00093 00094 double TripleRand::flat() { 00095 unsigned int ic ( integerCong ); 00096 unsigned int t ( tausworthe ); 00097 unsigned int h ( hurd ); 00098 return ( (t ^ ic ^ h) * twoToMinus_32() + // most significant part 00099 (h >> 11) * twoToMinus_53() + // fill in remaining bits 00100 nearlyTwoToMinus_54() // make sure non-zero 00101 ); 00102 } 00103 00104 void TripleRand::flatArray(const int size, double* vect) { 00105 for (int i = 0; i < size; ++i) { 00106 vect[i] = flat(); 00107 } 00108 } 00109 00110 void TripleRand::setSeed(long seed, int) { 00111 theSeed = seed; 00112 tausworthe = Tausworthe((unsigned int)seed + numEngines + 175321); 00113 integerCong = IntegerCong(69607 * tausworthe + 54329, numEngines); 00114 hurd = Hurd288Engine( 19781127 + integerCong ); 00115 } 00116 00117 void TripleRand::setSeeds(const long * seeds, int) { 00118 setSeed(seeds ? *seeds : 1234567, 0); 00119 theSeeds = seeds; 00120 } 00121 00122 void TripleRand::saveStatus(const char filename[]) const { 00123 std::ofstream outFile(filename, std::ios::out); 00124 if (!outFile.bad()) { 00125 outFile << "Uvec\n"; 00126 std::vector<unsigned long> v = put(); 00127 #ifdef TRACE_IO 00128 std::cout << "Result of v = put() is:\n"; 00129 #endif 00130 for (unsigned int i=0; i<v.size(); ++i) { 00131 outFile << v[i] << "\n"; 00132 #ifdef TRACE_IO 00133 std::cout << v[i] << " "; 00134 if (i%6==0) std::cout << "\n"; 00135 #endif 00136 } 00137 #ifdef TRACE_IO 00138 std::cout << "\n"; 00139 #endif 00140 } 00141 #ifdef REMOVED 00142 outFile << std::setprecision(20) << theSeed << " "; 00143 tausworthe.put ( outFile ); 00144 integerCong.put( outFile); 00145 outFile << ConstHurd() << std::endl; 00146 #endif 00147 } 00148 00149 void TripleRand::restoreStatus(const char filename[]) { 00150 std::ifstream inFile(filename, std::ios::in); 00151 if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) { 00152 std::cerr << " -- Engine state remains unchanged\n"; 00153 return; 00154 } 00155 if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) { 00156 std::vector<unsigned long> v; 00157 unsigned long xin; 00158 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) { 00159 inFile >> xin; 00160 #ifdef TRACE_IO 00161 std::cout << "ivec = " << ivec << " xin = " << xin << " "; 00162 if (ivec%3 == 0) std::cout << "\n"; 00163 #endif 00164 if (!inFile) { 00165 inFile.clear(std::ios::badbit | inFile.rdstate()); 00166 std::cerr << "\nTripleRand state (vector) description improper." 00167 << "\nrestoreStatus has failed." 00168 << "\nInput stream is probably mispositioned now." << std::endl; 00169 return; 00170 } 00171 v.push_back(xin); 00172 } 00173 getState(v); 00174 return; 00175 } 00176 00177 if (!inFile.bad()) { 00178 // inFile >> theSeed; removed -- encompased by possibleKeywordInput 00179 tausworthe.get ( inFile ); 00180 integerCong.get( inFile ); 00181 inFile >> Hurd(); 00182 } 00183 } 00184 00185 void TripleRand::showStatus() const { 00186 std::cout << std::setprecision(20) << std::endl; 00187 std::cout << "-------- TripleRand engine status ---------" 00188 << std::endl; 00189 std::cout << "Initial seed = " << theSeed << std::endl; 00190 std::cout << "Tausworthe generator = " << std::endl; 00191 tausworthe.put( std::cout ); 00192 std::cout << "IntegerCong generator = " << std::endl; 00193 integerCong.put( std::cout ); 00194 std::cout << "Hurd288Engine generator= " << std::endl << ConstHurd(); 00195 std::cout << std::endl << "-----------------------------------------" 00196 << std::endl; 00197 } 00198 00199 TripleRand::operator float() { 00200 return (float) 00201 ( ( integerCong ^ tausworthe ^ (unsigned int)hurd ) * twoToMinus_32() 00202 + nearlyTwoToMinus_54() ); 00203 // make sure non-zero! 00204 } 00205 00206 TripleRand::operator unsigned int() { 00207 return integerCong ^ tausworthe ^ (unsigned int)hurd; 00208 } 00209 00210 Hurd288Engine & TripleRand::Hurd() { return hurd; } 00211 00212 const Hurd288Engine & TripleRand::ConstHurd() const 00213 { return hurd; } 00214 00215 std::ostream & TripleRand::put (std::ostream & os ) const { 00216 char beginMarker[] = "TripleRand-begin"; 00217 os << beginMarker << "\nUvec\n"; 00218 std::vector<unsigned long> v = put(); 00219 for (unsigned int i=0; i<v.size(); ++i) { 00220 os << v[i] << "\n"; 00221 } 00222 return os; 00223 #ifdef REMOVED 00224 char endMarker[] = "TripleRand-end"; 00225 int pr=os.precision(20); 00226 os << " " << beginMarker << "\n"; 00227 os << theSeed << "\n"; 00228 tausworthe.put( os ); 00229 integerCong.put( os ); 00230 os << ConstHurd(); 00231 os << " " << endMarker << "\n"; 00232 os.precision(pr); 00233 return os; 00234 #endif 00235 } 00236 00237 std::vector<unsigned long> TripleRand::put () const { 00238 std::vector<unsigned long> v; 00239 v.push_back (engineIDulong<TripleRand>()); 00240 tausworthe.put(v); 00241 integerCong.put(v); 00242 std::vector<unsigned long> vHurd = hurd.put(); 00243 for (unsigned int i = 0; i < vHurd.size(); ++i) { 00244 v.push_back (vHurd[i]); 00245 } 00246 return v; 00247 } 00248 00249 std::istream & TripleRand::get (std::istream & is) { 00250 char beginMarker [MarkerLen]; 00251 is >> std::ws; 00252 is.width(MarkerLen); // causes the next read to the char* to be <= 00253 // that many bytes, INCLUDING A TERMINATION \0 00254 // (Stroustrup, section 21.3.2) 00255 is >> beginMarker; 00256 if (strcmp(beginMarker,"TripleRand-begin")) { 00257 is.clear(std::ios::badbit | is.rdstate()); 00258 std::cerr << "\nInput mispositioned or" 00259 << "\nTripleRand state description missing or" 00260 << "\nwrong engine type found." << std::endl; 00261 return is; 00262 } 00263 return getState(is); 00264 } 00265 00266 std::string TripleRand::beginTag ( ) { 00267 return "TripleRand-begin"; 00268 } 00269 00270 std::istream & TripleRand::getState (std::istream & is) { 00271 if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) { 00272 std::vector<unsigned long> v; 00273 unsigned long uu; 00274 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) { 00275 is >> uu; 00276 if (!is) { 00277 is.clear(std::ios::badbit | is.rdstate()); 00278 std::cerr << "\nTripleRand state (vector) description improper." 00279 << "\ngetState() has failed." 00280 << "\nInput stream is probably mispositioned now." << std::endl; 00281 return is; 00282 } 00283 v.push_back(uu); 00284 } 00285 getState(v); 00286 return (is); 00287 } 00288 00289 // is >> theSeed; Removed, encompassed by possibleKeywordInput() 00290 00291 char endMarker [MarkerLen]; 00292 tausworthe.get( is ); 00293 integerCong.get( is ); 00294 is >> Hurd(); 00295 is >> std::ws; 00296 is.width(MarkerLen); 00297 is >> endMarker; 00298 if (strcmp(endMarker,"TripleRand-end")) { 00299 is.clear(std::ios::badbit | is.rdstate()); 00300 std::cerr << "\nTripleRand state description incomplete." 00301 << "\nInput stream is probably mispositioned now." << std::endl; 00302 return is; 00303 } 00304 return is; 00305 } 00306 00307 bool TripleRand::get (const std::vector<unsigned long> & v) { 00308 if ((v[0] & 0xffffffffUL) != engineIDulong<TripleRand>()) { 00309 std::cerr << 00310 "\nTripleRand get:state vector has wrong ID word - state unchanged\n"; 00311 return false; 00312 } 00313 if (v.size() != VECTOR_STATE_SIZE) { 00314 std::cerr << "\nTripleRand get:state vector has wrong size: " 00315 << v.size() << " - state unchanged\n"; 00316 return false; 00317 } 00318 return getState(v); 00319 } 00320 00321 bool TripleRand::getState (const std::vector<unsigned long> & v) { 00322 std::vector<unsigned long>::const_iterator iv = v.begin()+1; 00323 if (!tausworthe.get(iv)) return false; 00324 if (!integerCong.get(iv)) return false; 00325 std::vector<unsigned long> vHurd; 00326 while (iv != v.end()) { 00327 vHurd.push_back(*iv++); 00328 } 00329 if (!hurd.get(vHurd)) { 00330 std::cerr << 00331 "\nTripleRand get from vector: problem getting the hurd sub-engine state\n"; 00332 return false; 00333 } 00334 return true; 00335 } 00336 00337 //******************************************************************** 00338 // Tausworthe 00339 //******************************************************************** 00340 00341 TripleRand::Tausworthe::Tausworthe() { 00342 words[0] = 1234567; 00343 for (wordIndex = 1; wordIndex < 4; ++wordIndex) { 00344 words[wordIndex] = 69607 * words[wordIndex-1] + 54329; 00345 } 00346 } 00347 00348 TripleRand::Tausworthe::Tausworthe(unsigned int seed) { 00349 words[0] = seed; 00350 for (wordIndex = 1; wordIndex < 4; ++wordIndex) { 00351 words[wordIndex] = 69607 * words[wordIndex-1] + 54329; 00352 } 00353 } 00354 00355 TripleRand::Tausworthe::operator unsigned int() { 00356 00357 // Mathematically: Consider a sequence of bits b[n]. Repeatedly form 00358 // b[0]' = b[127] ^ b[97]; b[n]' = b[n-1]. This sequence will have a very 00359 // long period (2**127-1 according to Tausworthe's work). 00360 00361 // The actual method used relies on the fact that the operations needed to 00362 // form bit 0 for up to 96 iterations never depend on the results of the 00363 // previous ones. So you can actually compute many bits at once. In fact 00364 // you can compute 32 at once -- despite 127 - 97 < 32 -- but 24 was used in 00365 // the method used in Canopy, where they wanted only single-precision float 00366 // randoms. I will do 32 here. 00367 00368 // When you do it this way, this looks disturbingly like the dread lagged XOR 00369 // Fibonacci. And indeed, it is a lagged Fibonacii, F(4,3, op) with the op 00370 // being the XOR of a combination of shifts of the two numbers. Although 00371 // Tausworthe asserted excellent properties, I would be scared to death. 00372 // However, the shifting and bit swapping really does randomize this in a 00373 // serious way. 00374 00375 // Statements have been made to the effect that shift register sequences fail 00376 // the parking lot test because they achieve randomness by multiple foldings, 00377 // and this produces a characteristic pattern. We observe that in this 00378 // specific algorithm, which has a fairly long lever arm, the foldings become 00379 // effectively random. This is evidenced by the fact that the generator 00380 // does pass the Diehard tests, including the parking lot test. 00381 00382 // To avoid shuffling of variables in memory, you either have to use circular 00383 // pointers (and those give you ifs, which are also costly) or compute at least 00384 // a few iterations at once. We do the latter. Although there is a possible 00385 // trade of room for more speed, by computing and saving 256 instead of 128 00386 // bits at once, I will stop at this level of optimization. 00387 00388 // To remind: Each (32-bit) step takes the XOR of bits [127-96] with bits 00389 // [95-64] and places it in bits [0-31]. But in the first step, we designate 00390 // word0 as bits [0-31], in the second step, word 1 (since the bits it holds 00391 // will no longer be needed), then word 2, then word 3. After this, the 00392 // stream contains 128 random bits which we will use as 4 valid 32-bit 00393 // random numbers. 00394 00395 // Thus at the start of the first step, word[0] contains the newest (used) 00396 // 32-bit random, and word[3] the oldest. After four steps, word[0] again 00397 // contains the newest (now unused) random, and word[3] the oldest. 00398 // Bit 0 of word[0] is logically the newest bit, and bit 31 of word[3] 00399 // the oldest. 00400 00401 if (wordIndex <= 0) { 00402 for (wordIndex = 0; wordIndex < 4; ++wordIndex) { 00403 words[wordIndex] = ( (words[(wordIndex+1) & 3] << 1 ) | 00404 (words[wordIndex] >> 31) ) 00405 ^ ( (words[(wordIndex+1) & 3] << 31) | 00406 (words[wordIndex] >> 1) ); 00407 } 00408 } 00409 return words[--wordIndex] & 0xffffffff; 00410 } 00411 00412 void TripleRand::Tausworthe::put( std::ostream & os ) const { 00413 char beginMarker[] = "Tausworthe-begin"; 00414 char endMarker[] = "Tausworthe-end"; 00415 00416 int pr=os.precision(20); 00417 os << " " << beginMarker << " "; 00418 os << std::setprecision(20); 00419 for (int i = 0; i < 4; ++i) { 00420 os << words[i] << " "; 00421 } 00422 os << wordIndex; 00423 os << " " << endMarker << " "; 00424 os << std::endl; 00425 os.precision(pr); 00426 } 00427 00428 void TripleRand::Tausworthe::put(std::vector<unsigned long> & v) const { 00429 for (int i = 0; i < 4; ++i) { 00430 v.push_back(static_cast<unsigned long>(words[i])); 00431 } 00432 v.push_back(static_cast<unsigned long>(wordIndex)); 00433 } 00434 00435 void TripleRand::Tausworthe::get( std::istream & is ) { 00436 char beginMarker [MarkerLen]; 00437 char endMarker [MarkerLen]; 00438 00439 is >> std::ws; 00440 is.width(MarkerLen); 00441 is >> beginMarker; 00442 if (strcmp(beginMarker,"Tausworthe-begin")) { 00443 is.clear(std::ios::badbit | is.rdstate()); 00444 std::cerr << "\nInput mispositioned or" 00445 << "\nTausworthe state description missing or" 00446 << "\nwrong engine type found." << std::endl; 00447 } 00448 for (int i = 0; i < 4; ++i) { 00449 is >> words[i]; 00450 } 00451 is >> wordIndex; 00452 is >> std::ws; 00453 is.width(MarkerLen); 00454 is >> endMarker; 00455 if (strcmp(endMarker,"Tausworthe-end")) { 00456 is.clear(std::ios::badbit | is.rdstate()); 00457 std::cerr << "\nTausworthe state description incomplete." 00458 << "\nInput stream is probably mispositioned now." << std::endl; 00459 } 00460 } 00461 00462 bool 00463 TripleRand::Tausworthe::get(std::vector<unsigned long>::const_iterator & iv){ 00464 for (int i = 0; i < 4; ++i) { 00465 words[i] = *iv++; 00466 } 00467 wordIndex = *iv++; 00468 return true; 00469 } 00470 00471 //******************************************************************** 00472 // IntegerCong 00473 //******************************************************************** 00474 00475 TripleRand::IntegerCong::IntegerCong() 00476 : state((unsigned int)3758656018U), 00477 multiplier(66565), 00478 addend(12341) 00479 { 00480 } 00481 00482 TripleRand::IntegerCong::IntegerCong(unsigned int seed, int streamNumber) 00483 : state(seed), 00484 multiplier(65536 + 1024 + 5 + (8 * 1017 * streamNumber)), 00485 addend(12341) 00486 { 00487 // As to the multiplier, the following comment was made: 00488 // We want our multipliers larger than 2^16, and equal to 00489 // 1 mod 4 (for max. period), but not equal to 1 mod 8 00490 // (for max. potency -- the smaller and higher dimension the 00491 // stripes, the better) 00492 00493 // All of these will have fairly long periods. Depending on the choice 00494 // of stream number, some of these will be quite decent when considered 00495 // as independent random engines, while others will be poor. Thus these 00496 // should not be used as stand-alone engines; but when combined with a 00497 // generator known to be good, they allow creation of millions of good 00498 // independent streams, without fear of two streams accidentally hitting 00499 // nearby places in the good random sequence. 00500 } 00501 00502 TripleRand::IntegerCong::operator unsigned int() { 00503 return state = (state * multiplier + addend) & 0xffffffff; 00504 } 00505 00506 void TripleRand::IntegerCong::put( std::ostream & os ) const { 00507 char beginMarker[] = "IntegerCong-begin"; 00508 char endMarker[] = "IntegerCong-end"; 00509 00510 int pr=os.precision(20); 00511 os << " " << beginMarker << " "; 00512 os << state << " " << multiplier << " " << addend; 00513 os << " " << endMarker << " "; 00514 os << std::endl; 00515 os.precision(pr); 00516 } 00517 00518 void TripleRand::IntegerCong::put(std::vector<unsigned long> & v) const { 00519 v.push_back(static_cast<unsigned long>(state)); 00520 v.push_back(static_cast<unsigned long>(multiplier)); 00521 v.push_back(static_cast<unsigned long>(addend)); 00522 } 00523 00524 void TripleRand::IntegerCong::get( std::istream & is ) { 00525 char beginMarker [MarkerLen]; 00526 char endMarker [MarkerLen]; 00527 00528 is >> std::ws; 00529 is.width(MarkerLen); 00530 is >> beginMarker; 00531 if (strcmp(beginMarker,"IntegerCong-begin")) { 00532 is.clear(std::ios::badbit | is.rdstate()); 00533 std::cerr << "\nInput mispositioned or" 00534 << "\nIntegerCong state description missing or" 00535 << "\nwrong engine type found." << std::endl; 00536 } 00537 is >> state >> multiplier >> addend; 00538 is >> std::ws; 00539 is.width(MarkerLen); 00540 is >> endMarker; 00541 if (strcmp(endMarker,"IntegerCong-end")) { 00542 is.clear(std::ios::badbit | is.rdstate()); 00543 std::cerr << "\nIntegerCong state description incomplete." 00544 << "\nInput stream is probably mispositioned now." << std::endl; 00545 } 00546 } 00547 00548 bool 00549 TripleRand::IntegerCong::get(std::vector<unsigned long>::const_iterator & iv) { 00550 state = *iv++; 00551 multiplier = *iv++; 00552 addend = *iv++; 00553 return true; 00554 } 00555 00556 } // namespace CLHEP