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

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