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