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

Hurd160Engine.cc
Go to the documentation of this file.
00001 // $Id: Hurd160Engine.cc,v 1.7 2010/07/20 18:07:17 garren Exp $
00002 // -*- C++ -*-
00003 //
00004 // -----------------------------------------------------------------------
00005 //                           HEP Random
00006 // Ken Smith      - Initial draft started:                    23rd Jul 1998
00007 //                - Added conversion operators:                6th Aug 1998
00008 // J. Marraffino  - Added some explicit casts to deal with
00009 //                  machines where sizeof(int) != sizeof(long)  22 Aug 1998
00010 // M. Fischler    - Modified use of the various exponents of 2
00011 //                  to avoid per-instance space overhead and
00012 //                  correct the rounding procedure              15 Sep 1998
00013 //                - Modified various methods to avoid any
00014 //                  possibility of predicting the next number
00015 //                  based on the last several:  We skip one 
00016 //                  32-bit word per cycle.                      15 Sep 1998
00017 //                - modify word[0] in two of the constructors
00018 //                  so that no sequence can ever accidentally
00019 //                  be produced by differnt seeds.              15 Sep 1998
00020 // J. Marraffino  - Remove dependence on hepStrings class       13 May 1999
00021 // M. Fischler    - Put endl at end of a save                   10 Apr 2001
00022 // M. Fischler    - In restore, checkFile for file not found    03 Dec 2004
00023 // M. Fischler    - Methods put, get for instance save/restore   12/8/04    
00024 // M. Fischler    - split get() into tag validation and 
00025 //                  getState() for anonymous restores           12/27/04    
00026 // M. Fischler    - put/get for vectors of ulongs               3/14/05
00027 // M. Fischler    - State-saving using only ints, for portability 4/12/05
00028 //                  
00029 // =======================================================================
00030 
00031 #include "CLHEP/Random/defs.h"
00032 #include "CLHEP/Random/Random.h"
00033 #include "CLHEP/Random/Hurd160Engine.h"
00034 #include "CLHEP/Random/engineIDulong.h"
00035 #include <string.h>     // for strcmp
00036 #include <cstdlib>      // for abs(int)
00037 
00038 using namespace std;
00039 
00040 namespace CLHEP {
00041 
00042 static const int MarkerLen = 64; // Enough room to hold a begin or end marker. 
00043 
00044 std::string Hurd160Engine::name() const {return "Hurd160Engine";}
00045 
00046 // Number of instances with automatic seed selection
00047 int Hurd160Engine::numEngines = 0;
00048 
00049 // Maximum index into the seed table
00050 int Hurd160Engine::maxIndex = 215;
00051 
00052 static inline unsigned int f160(unsigned int a, unsigned int b, unsigned int c)
00053 {
00054   return ( ((b<<2) & 0x7c) | ((a<<2) & ~0x7c) | (a>>30) ) ^ ( (c<<1)|(c>>31) );
00055 }
00056 
00057 Hurd160Engine::Hurd160Engine()
00058 : HepRandomEngine()
00059 {
00060   int cycle    = abs(int(numEngines/maxIndex));
00061   int curIndex = abs(int(numEngines%maxIndex));
00062   long mask = ((cycle & 0x007fffff) << 8);
00063   long seedlist[2];
00064   HepRandom::getTheTableSeeds( seedlist, curIndex );
00065   seedlist[0] ^= mask;
00066   seedlist[1] = 0;
00067   setSeeds(seedlist, 0);
00068   words[0] ^= 0x1324abcd;        // To make unique vs long or two unsigned
00069   if (words[0]==0) words[0] = 1; // ints in the constructor
00070   ++numEngines;
00071   for( int i=0; i < 100; ++i ) flat();            // warm-up just a bit
00072 }
00073 
00074 Hurd160Engine::Hurd160Engine( std::istream& is )
00075 : HepRandomEngine()
00076 {
00077     is >> *this;
00078 }
00079 
00080 Hurd160Engine::Hurd160Engine( long seed )
00081 : HepRandomEngine()
00082 {
00083   long seedlist[2]={seed,0};
00084   setSeeds(seedlist, 0);
00085   words[0] ^= 0xa5482134;        // To make unique vs default two unsigned
00086   if (words[0]==0) words[0] = 1; // ints in the constructor
00087   for( int i=0; i < 100; ++i ) flat();            // warm-up just a bit
00088 }
00089 
00090 Hurd160Engine::Hurd160Engine( int rowIndex, int colIndex )
00091 : HepRandomEngine()
00092 {
00093   int cycle = abs(int(rowIndex/maxIndex));
00094   int   row = abs(int(rowIndex%maxIndex));
00095   int   col = colIndex & 0x1;
00096   long mask = (( cycle & 0x000007ff ) << 20 );
00097   long seedlist[2];
00098   HepRandom::getTheTableSeeds( seedlist, row );
00099   // NOTE: is that really the seed wanted (PGC) ??
00100   seedlist[0] = (seedlist[col])^mask;
00101   seedlist[1]= 0;
00102   setSeeds(seedlist, 0);
00103   for( int i=0; i < 100; ++i ) flat();            // warm-up just a bit
00104 }
00105 
00106 Hurd160Engine::~Hurd160Engine() { }
00107 
00108 void Hurd160Engine::advance() {
00109      register unsigned int W0, W1, W2, W3, W4;
00110 
00111      W0 = words[0];
00112      W1 = words[1];
00113      W2 = words[2];
00114      W3 = words[3];
00115      W4 = words[4];
00116      W1 ^= W0; W0 = f160(W4, W3, W0);
00117      W2 ^= W1; W1 = f160(W0, W4, W1);
00118      W3 ^= W2; W2 = f160(W1, W0, W2);
00119      W4 ^= W3; W3 = f160(W2, W1, W3);
00120      W0 ^= W4; W4 = f160(W3, W2, W4);
00121      words[0] = W0 & 0xffffffff;
00122      words[1] = W1 & 0xffffffff;
00123      words[2] = W2 & 0xffffffff;
00124      words[3] = W3 & 0xffffffff;
00125      words[4] = W4 & 0xffffffff;
00126      wordIndex = 5;
00127 
00128 } // advance();
00129 
00130 double Hurd160Engine::flat() {
00131 
00132   if( wordIndex <= 2 ) {        // MF 9/15/98:
00133                                 // skip word 0 and use two words per flat
00134     advance();
00135   }
00136 
00137   // LG 6/30/2010
00138   // define the order of execution for --wordIndex
00139   double x = words[--wordIndex] * twoToMinus_32() ; // most significant part
00140   double y = (words[--wordIndex]>>11) * twoToMinus_53() + // fill in rest of bits
00141              nearlyTwoToMinus_54();        // make sure non-zero
00142   return  x + y ;
00143 }
00144 
00145 void Hurd160Engine::flatArray( const int size, double* vect ) {
00146     for (int i = 0; i < size; ++i) {
00147         vect[i] = flat();
00148     }
00149 }
00150 
00151 void Hurd160Engine::setSeed( long seed, int ) {
00152   words[0] = (unsigned int)seed;
00153   for (wordIndex = 1; wordIndex < 5; ++wordIndex) {
00154     words[wordIndex] = 69607 * words[wordIndex-1] + 54329;
00155   }
00156 }
00157 
00158 void Hurd160Engine::setSeeds( const long* seeds, int ) {
00159   setSeed( *seeds ? *seeds : 32767, 0 );
00160   theSeeds = seeds;
00161 }
00162      
00163 void Hurd160Engine::saveStatus( const char filename[] ) const {
00164   std::ofstream outFile(filename, std::ios::out);
00165   if( !outFile.bad() ) {
00166     outFile << "Uvec\n";
00167     std::vector<unsigned long> v = put();
00168                      #ifdef TRACE_IO
00169                          std::cout << "Result of v = put() is:\n"; 
00170                      #endif
00171     for (unsigned int i=0; i<v.size(); ++i) {
00172       outFile << v[i] << "\n";
00173                      #ifdef TRACE_IO
00174                            std::cout << v[i] << " ";
00175                            if (i%6==0) std::cout << "\n";
00176                      #endif
00177     }
00178                      #ifdef TRACE_IO
00179                          std::cout << "\n";
00180                      #endif
00181   }
00182 #ifdef REMOVED
00183     outFile << std::setprecision(20) << theSeed << " ";
00184     outFile << wordIndex << " ";
00185     for( int i = 0; i < 5; ++i ) {
00186       outFile << words[i] << " ";
00187     }
00188     outFile << std::endl;
00189 #endif
00190 }
00191 
00192 void Hurd160Engine::restoreStatus( const char filename[] ) {
00193   std::ifstream inFile(filename, std::ios::in);
00194   if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) { 
00195     std::cerr << "  -- Engine state remains unchanged\n";                 
00196     return;                                                               
00197   }                                                                       
00198   if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) {
00199     std::vector<unsigned long> v;
00200     unsigned long xin;
00201     for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
00202       inFile >> xin;
00203                #ifdef TRACE_IO
00204                std::cout << "ivec = " << ivec << "  xin = " << xin << "    ";
00205                if (ivec%3 == 0) std::cout << "\n"; 
00206                #endif
00207       if (!inFile) {
00208         inFile.clear(std::ios::badbit | inFile.rdstate());
00209         std::cerr << "\nHurd160Engine state (vector) description improper."
00210                << "\nrestoreStatus has failed."
00211                << "\nInput stream is probably mispositioned now." << std::endl;
00212         return;
00213       }
00214       v.push_back(xin);
00215     }
00216     getState(v);
00217     return;
00218   }
00219 
00220   if( !inFile.bad() ) {
00221 //     inFile >> theSeed;  removed -- encompased by possibleKeywordInput
00222     inFile >> wordIndex;
00223     for( int i = 0; i < 5; ++i ) {
00224         inFile >> words[i];
00225     }
00226   }
00227 }
00228 
00229 void Hurd160Engine::showStatus() const {
00230   int pr = std::cout.precision(20);
00231   std::cout << std::endl;
00232   std::cout << "----------- Hurd engine status ----------" << std::endl;
00233   std::cout << "Initial seed  = " << theSeed   << std::endl;
00234   std::cout << "Current index = " << wordIndex << std::endl;
00235   std::cout << "Current words = " << std::endl;
00236   for( int i = 0; i < 5 ; ++i ) {
00237     std::cout << "    " << words[i] << std::endl;
00238   }
00239   std::cout << "------------------------------------------" << std::endl;
00240   std::cout.precision(pr);
00241 }
00242 
00243 Hurd160Engine::operator float() {
00244   if( wordIndex <= 1 ) {        // MF 9/15/98:  skip word 0
00245     advance();
00246   }
00247   return words[--wordIndex ] * twoToMinus_32();
00248 }
00249 
00250 Hurd160Engine::operator unsigned int() {
00251   if( wordIndex <= 1 ) {        // MF 9/15/98:  skip word 0
00252     advance();
00253   }
00254   return words[--wordIndex];
00255 }
00256 
00257 std::ostream& Hurd160Engine::put(std::ostream& os) const {
00258   char beginMarker[] = "Hurd160Engine-begin";
00259   os << beginMarker << "\nUvec\n";
00260   std::vector<unsigned long> v = put();
00261   for (unsigned int i=0; i<v.size(); ++i) {
00262      os <<  v[i] <<  "\n";
00263   }
00264   return os;  
00265 #ifdef REMOVED 
00266   char endMarker[]   = "Hurd160Engine-end";
00267   int pr = os.precision(20);
00268   os << " " << beginMarker << " ";
00269   os << theSeed  << " ";
00270   os << wordIndex << " ";
00271   for (int i = 0; i < 5; ++i) {
00272     os << words[i]  << "\n";
00273   }
00274   os << endMarker   << "\n ";
00275   os.precision(pr);
00276   return os;
00277 #endif
00278 }
00279 
00280 std::vector<unsigned long> Hurd160Engine::put () const {
00281   std::vector<unsigned long> v;
00282   v.push_back (engineIDulong<Hurd160Engine>());
00283   v.push_back(static_cast<unsigned long>(wordIndex));
00284   for (int i = 0; i < 5; ++i) {
00285     v.push_back(static_cast<unsigned long>(words[i]));
00286   }
00287   return v;
00288 }
00289 
00290 
00291 std::istream& Hurd160Engine::get(std::istream& is) {
00292   char beginMarker [MarkerLen];
00293   is >> std::ws;
00294   is.width(MarkerLen);  // causes the next read to the char* to be <=
00295                         // that many bytes, INCLUDING A TERMINATION \0 
00296                         // (Stroustrup, section 21.3.2)
00297   is >> beginMarker;
00298   if (strcmp(beginMarker,"Hurd160Engine-begin")) {
00299     is.clear(std::ios::badbit | is.rdstate());
00300     std::cerr << "\nInput mispositioned or"
00301               << "\nHurd160Engine state description missing or"
00302               << "\nwrong engine type found." << std::endl;
00303     return is;
00304   }
00305   return getState(is);
00306 }
00307 
00308 std::string Hurd160Engine::beginTag ( )  { 
00309   return "Hurd160Engine-begin"; 
00310 }
00311 
00312 std::istream& Hurd160Engine::getState(std::istream& is) {
00313   if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) {
00314     std::vector<unsigned long> v;
00315     unsigned long uu;
00316     for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
00317       is >> uu;
00318       if (!is) {
00319         is.clear(std::ios::badbit | is.rdstate());
00320         std::cerr << "\nHurd160Engine state (vector) description improper."
00321                 << "\ngetState() has failed."
00322                << "\nInput stream is probably mispositioned now." << std::endl;
00323         return is;
00324       }
00325       v.push_back(uu);
00326     }
00327     getState(v);
00328     return (is);
00329   }
00330 
00331 //  is >> theSeed;  Removed, encompassed by possibleKeywordInput()
00332 
00333   char endMarker   [MarkerLen];
00334   is >> wordIndex;
00335   for (int i = 0; i < 5; ++i) {
00336     is >> words[i];
00337   }
00338   is >> std::ws;
00339   is.width(MarkerLen);
00340   is >> endMarker;
00341   if (strcmp(endMarker,"Hurd160Engine-end")) {
00342     is.clear(std::ios::badbit | is.rdstate());
00343     std::cerr << "\nHurd160Engine state description incomplete."
00344               << "\nInput stream is probably mispositioned now." << std::endl;
00345     return is;
00346   }
00347   return is;
00348 }
00349 
00350 
00351 bool Hurd160Engine::get (const std::vector<unsigned long> & v) {
00352   if ((v[0] & 0xffffffffUL) != engineIDulong<Hurd160Engine>()) {
00353     std::cerr << 
00354         "\nHurd160Engine get:state vector has wrong ID word - state unchanged\n";
00355     return false;
00356   }
00357   return getState(v);
00358 }
00359 
00360 bool Hurd160Engine::getState (const std::vector<unsigned long> & v) {
00361   if (v.size() != VECTOR_STATE_SIZE ) {
00362     std::cerr << 
00363         "\nHurd160Engine get:state vector has wrong length - state unchanged\n";
00364     return false;
00365   }
00366   wordIndex = v[1];
00367   for (int i = 0; i < 5; ++i) {
00368     words[i] = v[i+2];
00369   }
00370   return true;
00371 }
00372 
00373 
00374 }  // namespace CLHEP