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