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