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