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