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