CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: Hurd288Engine.cc,v 1.4.4.2.2.2 2010/06/30 17:11:00 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> 00042 #include <cmath> // for ldexp() 00043 #include <stdlib.h> // for abs(int) 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 double Hurd288Engine::twoToMinus_32; 00052 double Hurd288Engine::twoToMinus_53; 00053 double Hurd288Engine::nearlyTwoToMinus_54; 00054 00055 std::string Hurd288Engine::name() const {return "Hurd288Engine";} 00056 00057 void Hurd288Engine::powersOfTwo() { 00058 twoToMinus_32 = ldexp (1.0, -32); 00059 twoToMinus_53 = ldexp (1.0, -53); 00060 nearlyTwoToMinus_54 = ldexp (1.0, -54) - ldexp (1.0, -100); 00061 } 00062 00063 // Number of instances with automatic seed selection 00064 int Hurd288Engine::numEngines = 0; 00065 00066 // Maximum index into the seed table 00067 int Hurd288Engine::maxIndex = 215; 00068 00069 static inline unsigned int f288(unsigned int a, unsigned int b, unsigned int c) 00070 { 00071 return ( ((b<<2) & 0x7ffc) | ((a<<2) & ~0x7ffc) | (a>>30) ) ^ 00072 ( (c<<1)|(c>>31) ); 00073 } 00074 00075 Hurd288Engine::Hurd288Engine() { 00076 powersOfTwo(); 00077 int cycle = abs(int(numEngines/maxIndex)); 00078 int curIndex = abs(int(numEngines%maxIndex)); 00079 long mask = ((cycle & 0x007fffff) << 8); 00080 long seedlist[2]; 00081 HepRandom::getTheTableSeeds( seedlist, curIndex ); 00082 seedlist[0] ^= mask; 00083 seedlist[1] = 0; 00084 setSeeds(seedlist, 0); 00085 words[0] ^= 0x1324abcd; // To make unique vs long or two unsigned 00086 if (words[0]==0) words[0] = 1; // ints in the constructor 00087 ++numEngines; 00088 for( int i=0; i < 100; ++i ) flat(); // warm up just a bit 00089 } 00090 00091 Hurd288Engine::Hurd288Engine( std::istream& is ) 00092 { 00093 is >> *this; 00094 } 00095 00096 Hurd288Engine::Hurd288Engine( long seed ) { 00097 powersOfTwo(); 00098 long seedlist[2]={seed,0}; 00099 setSeeds(seedlist, 0); 00100 words[0] ^= 0xa5482134; // To make unique vs default two unsigned 00101 if (words[0]==0) words[0] = 1; // ints in the constructor 00102 for( int i=0; i < 100; ++i ) flat(); // warm up just a bit 00103 } 00104 00105 Hurd288Engine::Hurd288Engine( int rowIndex, int colIndex ) { 00106 powersOfTwo(); 00107 int cycle = abs(int(rowIndex/maxIndex)); 00108 int row = abs(int(rowIndex%maxIndex)); 00109 int col = colIndex & 0x1; 00110 long mask = (( cycle & 0x000007ff ) << 20 ); 00111 long seedlist[2]; 00112 HepRandom::getTheTableSeeds( seedlist, row ); 00113 seedlist[0] = (seedlist[col])^mask; 00114 seedlist[1]= 0; 00115 setSeeds(seedlist, 0); 00116 for( int i=0; i < 100; ++i ) flat(); // warm up just a bit 00117 } 00118 00119 Hurd288Engine::~Hurd288Engine() { } 00120 00121 Hurd288Engine::Hurd288Engine( const Hurd288Engine & p ) 00122 { 00123 *this = p; 00124 } 00125 00126 Hurd288Engine & Hurd288Engine::operator=( const Hurd288Engine & p ) { 00127 if( this != &p ) { 00128 wordIndex = p.wordIndex; 00129 for( int i=0; i < 9; ++i ) { 00130 words[i] = p.words[i]; 00131 } 00132 } 00133 return *this; 00134 } 00135 00136 void Hurd288Engine::advance() { 00137 00138 register unsigned int W0, W1, W2, W3, W4, W5, W6, W7, W8; 00139 00140 W0 = words[0]; 00141 W1 = words[1]; 00142 W2 = words[2]; 00143 W3 = words[3]; 00144 W4 = words[4]; 00145 W5 = words[5]; 00146 W6 = words[6]; 00147 W7 = words[7]; 00148 W8 = words[8]; 00149 W1 ^= W0; W0 = f288(W2, W3, W0); 00150 W2 ^= W1; W1 = f288(W3, W4, W1); 00151 W3 ^= W2; W2 = f288(W4, W5, W2); 00152 W4 ^= W3; W3 = f288(W5, W6, W3); 00153 W5 ^= W4; W4 = f288(W6, W7, W4); 00154 W6 ^= W5; W5 = f288(W7, W8, W5); 00155 W7 ^= W6; W6 = f288(W8, W0, W6); 00156 W8 ^= W7; W7 = f288(W0, W1, W7); 00157 W0 ^= W8; W8 = f288(W1, W2, W8); 00158 words[0] = W0 & 0xffffffff; 00159 words[1] = W1 & 0xffffffff; 00160 words[2] = W2 & 0xffffffff; 00161 words[3] = W3 & 0xffffffff; 00162 words[4] = W4 & 0xffffffff; 00163 words[5] = W5 & 0xffffffff; 00164 words[6] = W6 & 0xffffffff; 00165 words[7] = W7 & 0xffffffff; 00166 words[8] = W8 & 0xffffffff; 00167 wordIndex = 9; 00168 00169 } // advance() 00170 00171 00172 double Hurd288Engine::flat() { 00173 00174 if( wordIndex <= 2 ) { // MF 9/15/98: 00175 // skip word 0 and use two words per flat 00176 advance(); 00177 } 00178 00179 // LG 6/30/2010 00180 // define the order of execution for --wordIndex 00181 double x = words[--wordIndex] * twoToMinus_32 ; // most significant part 00182 double y = (words[--wordIndex]>>11) * twoToMinus_53 + // fill in rest of bits 00183 nearlyTwoToMinus_54; // make sure non-zero 00184 return x + y; 00185 } 00186 00187 void Hurd288Engine::flatArray( const int size, double* vect ) { 00188 for (int i = 0; i < size; ++i) { 00189 vect[i] = flat(); 00190 } 00191 } 00192 00193 void Hurd288Engine::setSeed( long seed, int ) { 00194 words[0] = (unsigned int)seed; 00195 for (wordIndex = 1; wordIndex < 9; ++wordIndex) { 00196 words[wordIndex] = 69607 * words[wordIndex-1] + 54329; 00197 } 00198 } 00199 00200 void Hurd288Engine::setSeeds( const long* seeds, int ) { 00201 setSeed( *seeds ? *seeds : 32767, 0 ); 00202 theSeeds = seeds; 00203 } 00204 00205 void Hurd288Engine::saveStatus( const char filename[] ) const { 00206 std::ofstream outFile(filename, std::ios::out); 00207 if( !outFile.bad() ) { 00208 outFile << "Uvec\n"; 00209 std::vector<unsigned long> v = put(); 00210 #ifdef TRACE_IO 00211 std::cout << "Result of v = put() is:\n"; 00212 #endif 00213 for (unsigned int i=0; i<v.size(); ++i) { 00214 outFile << v[i] << "\n"; 00215 #ifdef TRACE_IO 00216 std::cout << v[i] << " "; 00217 if (i%6==0) std::cout << "\n"; 00218 #endif 00219 } 00220 #ifdef TRACE_IO 00221 std::cout << "\n"; 00222 #endif 00223 } 00224 #ifdef REMOVED 00225 outFile << std::setprecision(20) << theSeed << " "; 00226 outFile << wordIndex << " "; 00227 for( int i = 0; i < 9; ++i ) { 00228 outFile << words[i] << " "; 00229 } 00230 outFile << std::endl; 00231 #endif 00232 } 00233 00234 void Hurd288Engine::restoreStatus( const char filename[] ) { 00235 std::ifstream inFile(filename, std::ios::in); 00236 if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) { 00237 std::cerr << " -- Engine state remains unchanged\n"; 00238 return; 00239 } 00240 if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) { 00241 std::vector<unsigned long> v; 00242 unsigned long xin; 00243 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) { 00244 inFile >> xin; 00245 #ifdef TRACE_IO 00246 std::cout << "ivec = " << ivec << " xin = " << xin << " "; 00247 if (ivec%3 == 0) std::cout << "\n"; 00248 #endif 00249 if (!inFile) { 00250 inFile.clear(std::ios::badbit | inFile.rdstate()); 00251 std::cerr << "\nHurd288Engine state (vector) description improper." 00252 << "\nrestoreStatus has failed." 00253 << "\nInput stream is probably mispositioned now." << std::endl; 00254 return; 00255 } 00256 v.push_back(xin); 00257 } 00258 getState(v); 00259 return; 00260 } 00261 00262 if( !inFile.bad() ) { 00263 // inFile >> theSeed; removed -- encompased by possibleKeywordInput 00264 inFile >> wordIndex; 00265 for( int i = 0; i < 9; ++i ) { 00266 inFile >> words[i]; 00267 } 00268 } 00269 } 00270 00271 void Hurd288Engine::showStatus() const { 00272 std::cout << std::setprecision(20) << std::endl; 00273 std::cout << "----------- Hurd2 engine status ----------" << std::endl; 00274 std::cout << "Initial seed = " << theSeed << std::endl; 00275 std::cout << "Current index = " << wordIndex << std::endl; 00276 std::cout << "Current words = " << std::endl; 00277 for( int i = 0; i < 9 ; ++i ) { 00278 std::cout << " " << words[i] << std::endl; 00279 } 00280 std::cout << "-------------------------------------------" << std::endl; 00281 } 00282 00283 Hurd288Engine::operator float() { 00284 if( wordIndex <= 1 ) { // MF 9/15/98: skip word 0 00285 advance(); 00286 } 00287 return words[--wordIndex ] * twoToMinus_32; 00288 } 00289 00290 Hurd288Engine::operator unsigned int() { 00291 if( wordIndex <= 1 ) { // MF 9/15/98: skip word 0 00292 advance(); 00293 } 00294 return words[--wordIndex]; 00295 } 00296 00297 std::ostream& Hurd288Engine::put(std::ostream& os) const { 00298 char beginMarker[] = "Hurd288Engine-begin"; 00299 os << beginMarker << "\nUvec\n"; 00300 std::vector<unsigned long> v = put(); 00301 for (unsigned int i=0; i<v.size(); ++i) { 00302 os << v[i] << "\n"; 00303 } 00304 return os; 00305 #ifdef REMOVED 00306 char endMarker[] = "Hurd288Engine-end"; 00307 int pr = os.precision(20); 00308 os << " " << beginMarker << " "; 00309 os << theSeed << " "; 00310 os << wordIndex << " "; 00311 for (int i = 0; i < 9; ++i) { 00312 os << words[i] << "\n"; 00313 } 00314 os << endMarker << "\n "; 00315 os.precision(pr); 00316 return os; 00317 #endif 00318 } 00319 00320 std::vector<unsigned long> Hurd288Engine::put () const { 00321 std::vector<unsigned long> v; 00322 v.push_back (engineIDulong<Hurd288Engine>()); 00323 v.push_back(static_cast<unsigned long>(wordIndex)); 00324 for (int i = 0; i < 9; ++i) { 00325 v.push_back(static_cast<unsigned long>(words[i])); 00326 } 00327 return v; 00328 } 00329 00330 00331 std::istream& Hurd288Engine::get(std::istream& is) { 00332 char beginMarker [MarkerLen]; 00333 is >> std::ws; 00334 is.width(MarkerLen); // causes the next read to the char* to be <= 00335 // that many bytes, INCLUDING A TERMINATION \0 00336 // (Stroustrup, section 21.3.2) 00337 is >> beginMarker; 00338 if (strcmp(beginMarker,"Hurd288Engine-begin")) { 00339 is.clear(std::ios::badbit | is.rdstate()); 00340 std::cerr << "\nInput mispositioned or" 00341 << "\nHurd288Engine state description missing or" 00342 << "\nwrong engine type found." << std::endl; 00343 return is; 00344 } 00345 return getState(is); 00346 } 00347 00348 std::string Hurd288Engine::beginTag ( ) { 00349 return "Hurd288Engine-begin"; 00350 } 00351 00352 std::istream& Hurd288Engine::getState(std::istream& is) { 00353 if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) { 00354 std::vector<unsigned long> v; 00355 unsigned long uu; 00356 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) { 00357 is >> uu; 00358 if (!is) { 00359 is.clear(std::ios::badbit | is.rdstate()); 00360 std::cerr << "\nHurd288Engine state (vector) description improper." 00361 << "\ngetState() has failed." 00362 << "\nInput stream is probably mispositioned now." << std::endl; 00363 return is; 00364 } 00365 v.push_back(uu); 00366 } 00367 getState(v); 00368 return (is); 00369 } 00370 00371 // is >> theSeed; Removed, encompassed by possibleKeywordInput() 00372 00373 char endMarker [MarkerLen]; 00374 is >> wordIndex; 00375 for (int i = 0; i < 9; ++i) { 00376 is >> words[i]; 00377 } 00378 is >> std::ws; 00379 is.width(MarkerLen); 00380 is >> endMarker; 00381 if (strcmp(endMarker,"Hurd288Engine-end")) { 00382 is.clear(std::ios::badbit | is.rdstate()); 00383 std::cerr << "\nHurd288Engine state description incomplete." 00384 << "\nInput stream is probably mispositioned now." << std::endl; 00385 return is; 00386 } 00387 return is; 00388 } 00389 00390 bool Hurd288Engine::get (const std::vector<unsigned long> & v) { 00391 if ((v[0] & 0xffffffffUL) != engineIDulong<Hurd288Engine>()) { 00392 std::cerr << 00393 "\nHurd288Engine get:state vector has wrong ID word - state unchanged\n"; 00394 std::cerr << "The correct ID would be " << engineIDulong<Hurd288Engine>() 00395 << "; the actual ID is " << v[0] << "\n"; 00396 return false; 00397 } 00398 return getState(v); 00399 } 00400 00401 bool Hurd288Engine::getState (const std::vector<unsigned long> & v) { 00402 if (v.size() != VECTOR_STATE_SIZE ) { 00403 std::cerr << 00404 "\nHurd288Engine get:state vector has wrong length - state unchanged\n"; 00405 return false; 00406 } 00407 wordIndex = v[1]; 00408 for (int i = 0; i < 9; ++i) { 00409 words[i] = v[i+2]; 00410 } 00411 return true; 00412 } 00413 00414 } // namespace CLHEP