CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
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