CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

DualRand.cc
Go to the documentation of this file.
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