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