CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: RanecuEngine.cc,v 1.4.4.2.2.8 2010/06/23 20:49:50 garren Exp $ 00002 // -*- C++ -*- 00003 // 00004 // ----------------------------------------------------------------------- 00005 // HEP Random 00006 // --- RanecuEngine --- 00007 // class implementation file 00008 // ----------------------------------------------------------------------- 00009 // This file is part of Geant4 (simulation toolkit for HEP). 00010 // 00011 // RANECU Random Engine - algorithm originally written in FORTRAN77 00012 // as part of the MATHLIB HEP library. 00013 00014 // ======================================================================= 00015 // Gabriele Cosmo - Created - 2nd February 1996 00016 // - Minor corrections: 31st October 1996 00017 // - Added methods for engine status: 19th November 1996 00018 // - Added abs for setting seed index: 11th July 1997 00019 // - Modified setSeeds() to handle default index: 16th Oct 1997 00020 // - setSeed() now resets the engine status to the original 00021 // values in the static table of HepRandom: 19th Mar 1998 00022 // J.Marraffino - Added stream operators and related constructor. 00023 // Added automatic seed selection from seed table and 00024 // engine counter: 16th Feb 1998 00025 // Ken Smith - Added conversion operators: 6th Aug 1998 00026 // J. Marraffino - Remove dependence on hepString class 13 May 1999 00027 // M. Fischler - Add endl to the end of saveStatus 10 Apr 2001 00028 // M. Fischler - In restore, checkFile for file not found 03 Dec 2004 00029 // M. Fischler - Methods for distrib. instance 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 // M. Fischler - Modify ctor and setSeed to utilize all info provided 00035 // and avoid coincidence of same state from different 00036 // seeds 6/22/10 00037 // 00038 // ======================================================================= 00039 00040 #include "CLHEP/Random/defs.h" 00041 #include "CLHEP/Random/Random.h" 00042 #include "CLHEP/Random/RanecuEngine.h" 00043 #include "CLHEP/Random/engineIDulong.h" 00044 #include <string.h> 00045 #include <cmath> 00046 #include <stdlib.h> 00047 00048 namespace CLHEP { 00049 00050 static const int MarkerLen = 64; // Enough room to hold a begin or end marker. 00051 00052 std::string RanecuEngine::name() const {return "RanecuEngine";} 00053 00054 void RanecuEngine::further_randomize (int seq, int col, int index, int modulus) 00055 { 00056 table[seq][col] -= (index&0x3FFFFFFF); 00057 while (table[seq][col] <= 0) table[seq][col] += (modulus-1); 00058 } // mf 6/22/10 00059 00060 // Number of instances with automatic seed selection 00061 int RanecuEngine::numEngines = 0; 00062 00063 RanecuEngine::RanecuEngine() 00064 : ecuyer_a(40014),ecuyer_b(53668),ecuyer_c(12211), 00065 ecuyer_d(40692),ecuyer_e(52774),ecuyer_f(3791), 00066 shift1(2147483563),shift2(2147483399), 00067 prec(4.6566128E-10 ) 00068 { 00069 int cycle = abs(int(numEngines/maxSeq)); 00070 seq = abs(int(numEngines%maxSeq)); 00071 numEngines += 1; 00072 theSeed = seq; 00073 long mask = ((cycle & 0x007fffff) << 8); 00074 for (int i=0; i<2; ++i) { 00075 for (int j=0; j<maxSeq; ++j) { 00076 HepRandom::getTheTableSeeds(table[j],j); 00077 table[j][i] ^= mask; 00078 } 00079 } 00080 theSeeds = &table[seq][0]; 00081 } 00082 00083 RanecuEngine::RanecuEngine(int index) 00084 : ecuyer_a(40014),ecuyer_b(53668),ecuyer_c(12211), 00085 ecuyer_d(40692),ecuyer_e(52774),ecuyer_f(3791), 00086 shift1(2147483563),shift2(2147483399), 00087 prec(4.6566128E-10 ) 00088 { 00089 int cycle = abs(int(index/maxSeq)); 00090 seq = abs(int(index%maxSeq)); 00091 theSeed = seq; 00092 long mask = ((cycle & 0x000007ff) << 20); 00093 for (int j=0; j<maxSeq; ++j) { 00094 HepRandom::getTheTableSeeds(table[j],j); 00095 table[j][0] ^= mask; 00096 table[j][1] ^= mask; 00097 } 00098 theSeeds = &table[seq][0]; 00099 further_randomize (seq, 0, index, shift1); // mf 6/22/10 00100 } 00101 00102 RanecuEngine::RanecuEngine(std::istream& is) 00103 : ecuyer_a(40014),ecuyer_b(53668),ecuyer_c(12211), 00104 ecuyer_d(40692),ecuyer_e(52774),ecuyer_f(3791), 00105 shift1(2147483563),shift2(2147483399), 00106 prec(4.6566128E-10 ) 00107 { 00108 is >> *this; 00109 } 00110 00111 RanecuEngine::~RanecuEngine() {} 00112 00113 RanecuEngine::RanecuEngine(const RanecuEngine &p) 00114 : ecuyer_a(40014),ecuyer_b(53668),ecuyer_c(12211), 00115 ecuyer_d(40692),ecuyer_e(52774),ecuyer_f(3791), 00116 shift1(2147483563),shift2(2147483399), 00117 prec(4.6566128E-10 ) 00118 { 00119 if ((this != &p) && (&p)) { 00120 theSeed = p.getSeed(); 00121 seq = p.seq; 00122 for (int i=0; i<2; ++i) 00123 for (int j=0; j<maxSeq; ++j) 00124 table[j][i] = p.table[j][i]; 00125 theSeeds = &table[seq][0]; 00126 } 00127 } 00128 00129 RanecuEngine & RanecuEngine::operator = (const RanecuEngine &p) 00130 { 00131 if ((this != &p) && (&p)) { 00132 theSeed = p.getSeed(); 00133 seq = p.seq; 00134 for (int i=0; i<2; ++i) 00135 for (int j=0; j<maxSeq; ++j) 00136 table[j][i] = p.table[j][i]; 00137 theSeeds = &table[seq][0]; 00138 } 00139 return *this; 00140 } 00141 00142 void RanecuEngine::setSeed(long index, int dum) 00143 { 00144 seq = abs(int(index%maxSeq)); 00145 theSeed = seq; 00146 HepRandom::getTheTableSeeds(table[seq],seq); 00147 theSeeds = &table[seq][0]; 00148 further_randomize (seq, 0, index, shift1); // mf 6/22/10 00149 further_randomize (seq, 1, dum, shift2); // mf 6/22/10 00150 } 00151 00152 void RanecuEngine::setSeeds(const long* seeds, int pos) 00153 { 00154 if (pos != -1) { 00155 seq = abs(int(pos%maxSeq)); 00156 theSeed = seq; 00157 } 00158 // only positive seeds are allowed 00159 table[seq][0] = abs(seeds[0])%shift1; 00160 table[seq][1] = abs(seeds[1])%shift2; 00161 theSeeds = &table[seq][0]; 00162 } 00163 00164 void RanecuEngine::setIndex(long index) 00165 { 00166 seq = abs(int(index%maxSeq)); 00167 theSeed = seq; 00168 theSeeds = &table[seq][0]; 00169 } 00170 00171 void RanecuEngine::saveStatus( const char filename[] ) const 00172 { 00173 std::ofstream outFile( filename, std::ios::out ) ; 00174 00175 if (!outFile.bad()) { 00176 outFile << "Uvec\n"; 00177 std::vector<unsigned long> v = put(); 00178 #ifdef TRACE_IO 00179 std::cout << "Result of v = put() is:\n"; 00180 #endif 00181 for (unsigned int i=0; i<v.size(); ++i) { 00182 outFile << v[i] << "\n"; 00183 #ifdef TRACE_IO 00184 std::cout << v[i] << " "; 00185 if (i%6==0) std::cout << "\n"; 00186 #endif 00187 } 00188 #ifdef TRACE_IO 00189 std::cout << "\n"; 00190 #endif 00191 } 00192 #ifdef REMOVED 00193 if (!outFile.bad()) { 00194 outFile << theSeed << std::endl; 00195 for (int i=0; i<2; ++i) 00196 outFile << table[theSeed][i] << " "; 00197 outFile << std::endl; 00198 } 00199 #endif 00200 } 00201 00202 void RanecuEngine::restoreStatus( const char filename[] ) 00203 { 00204 std::ifstream inFile( filename, std::ios::in); 00205 if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) { 00206 std::cerr << " -- Engine state remains unchanged\n"; 00207 return; 00208 } 00209 if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) { 00210 std::vector<unsigned long> v; 00211 unsigned long xin; 00212 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) { 00213 inFile >> xin; 00214 #ifdef TRACE_IO 00215 std::cout << "ivec = " << ivec << " xin = " << xin << " "; 00216 if (ivec%3 == 0) std::cout << "\n"; 00217 #endif 00218 if (!inFile) { 00219 inFile.clear(std::ios::badbit | inFile.rdstate()); 00220 std::cerr << "\nJamesRandom state (vector) description improper." 00221 << "\nrestoreStatus has failed." 00222 << "\nInput stream is probably mispositioned now." << std::endl; 00223 return; 00224 } 00225 v.push_back(xin); 00226 } 00227 getState(v); 00228 return; 00229 } 00230 00231 if (!inFile.bad() && !inFile.eof()) { 00232 // inFile >> theSeed; removed -- encompased by possibleKeywordInput 00233 for (int i=0; i<2; ++i) 00234 inFile >> table[theSeed][i]; 00235 seq = int(theSeed); 00236 } 00237 } 00238 00239 void RanecuEngine::showStatus() const 00240 { 00241 std::cout << std::endl; 00242 std::cout << "--------- Ranecu engine status ---------" << std::endl; 00243 std::cout << " Initial seed (index) = " << theSeed << std::endl; 00244 std::cout << " Current couple of seeds = " 00245 << table[theSeed][0] << ", " 00246 << table[theSeed][1] << std::endl; 00247 std::cout << "----------------------------------------" << std::endl; 00248 } 00249 00250 double RanecuEngine::flat() 00251 { 00252 const int index = seq; 00253 long seed1 = table[index][0]; 00254 long seed2 = table[index][1]; 00255 00256 int k1 = (int)(seed1/ecuyer_b); 00257 int k2 = (int)(seed2/ecuyer_e); 00258 00259 seed1 = ecuyer_a*(seed1-k1*ecuyer_b)-k1*ecuyer_c; 00260 if (seed1 < 0) seed1 += shift1; 00261 seed2 = ecuyer_d*(seed2-k2*ecuyer_e)-k2*ecuyer_f; 00262 if (seed2 < 0) seed2 += shift2; 00263 00264 table[index][0] = seed1; 00265 table[index][1] = seed2; 00266 00267 long diff = seed1-seed2; 00268 00269 if (diff <= 0) diff += (shift1-1); 00270 return (double)(diff*prec); 00271 } 00272 00273 void RanecuEngine::flatArray(const int size, double* vect) 00274 { 00275 const int index = seq; 00276 long seed1 = table[index][0]; 00277 long seed2 = table[index][1]; 00278 int k1, k2; 00279 register int i; 00280 00281 for (i=0; i<size; ++i) 00282 { 00283 k1 = (int)(seed1/ecuyer_b); 00284 k2 = (int)(seed2/ecuyer_e); 00285 00286 seed1 = ecuyer_a*(seed1-k1*ecuyer_b)-k1*ecuyer_c; 00287 if (seed1 < 0) seed1 += shift1; 00288 seed2 = ecuyer_d*(seed2-k2*ecuyer_e)-k2*ecuyer_f; 00289 if (seed2 < 0) seed2 += shift2; 00290 00291 long diff = seed1-seed2; 00292 if (diff <= 0) diff += (shift1-1); 00293 00294 vect[i] = (double)(diff*prec); 00295 } 00296 table[index][0] = seed1; 00297 table[index][1] = seed2; 00298 } 00299 00300 RanecuEngine::operator unsigned int() { 00301 const int index = seq; 00302 long seed1 = table[index][0]; 00303 long seed2 = table[index][1]; 00304 00305 int k1 = (int)(seed1/ecuyer_b); 00306 int k2 = (int)(seed2/ecuyer_e); 00307 00308 seed1 = ecuyer_a*(seed1-k1*ecuyer_b)-k1*ecuyer_c; 00309 if (seed1 < 0) seed1 += shift1; 00310 seed2 = ecuyer_d*(seed2-k2*ecuyer_e)-k2*ecuyer_f; 00311 if (seed2 < 0) seed2 += shift2; 00312 00313 table[index][0] = seed1; 00314 table[index][1] = seed2; 00315 long diff = seed1-seed2; 00316 if( diff <= 0 ) diff += (shift1-1); 00317 00318 return ((diff << 1) | (seed1&1))& 0xffffffff; 00319 } 00320 00321 std::ostream & RanecuEngine::put( std::ostream& os ) const 00322 { 00323 char beginMarker[] = "RanecuEngine-begin"; 00324 os << beginMarker << "\nUvec\n"; 00325 std::vector<unsigned long> v = put(); 00326 for (unsigned int i=0; i<v.size(); ++i) { 00327 os << v[i] << "\n"; 00328 } 00329 return os; 00330 #ifdef REMOVED 00331 char endMarker[] = "RanecuEngine-end"; 00332 os << " " << beginMarker << "\n"; 00333 os << theSeed << " "; 00334 for (int i=0; i<2; ++i) { 00335 os << table[theSeed][i] << "\n"; 00336 } 00337 os << endMarker << "\n"; 00338 return os; 00339 #endif 00340 } 00341 00342 std::vector<unsigned long> RanecuEngine::put () const { 00343 std::vector<unsigned long> v; 00344 v.push_back (engineIDulong<RanecuEngine>()); 00345 v.push_back(static_cast<unsigned long>(theSeed)); 00346 v.push_back(static_cast<unsigned long>(table[theSeed][0])); 00347 v.push_back(static_cast<unsigned long>(table[theSeed][1])); 00348 return v; 00349 } 00350 00351 std::istream & RanecuEngine::get ( std::istream& is ) 00352 { 00353 char beginMarker [MarkerLen]; 00354 00355 is >> std::ws; 00356 is.width(MarkerLen); // causes the next read to the char* to be <= 00357 // that many bytes, INCLUDING A TERMINATION \0 00358 // (Stroustrup, section 21.3.2) 00359 is >> beginMarker; 00360 if (strcmp(beginMarker,"RanecuEngine-begin")) { 00361 is.clear(std::ios::badbit | is.rdstate()); 00362 std::cerr << "\nInput stream mispositioned or" 00363 << "\nRanecuEngine state description missing or" 00364 << "\nwrong engine type found." << std::endl; 00365 return is; 00366 } 00367 return getState(is); 00368 } 00369 00370 std::string RanecuEngine::beginTag ( ) { 00371 return "RanecuEngine-begin"; 00372 } 00373 00374 std::istream & RanecuEngine::getState ( std::istream& is ) 00375 { 00376 if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) { 00377 std::vector<unsigned long> v; 00378 unsigned long uu; 00379 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) { 00380 is >> uu; 00381 if (!is) { 00382 is.clear(std::ios::badbit | is.rdstate()); 00383 std::cerr << "\nRanecuEngine state (vector) description improper." 00384 << "\ngetState() has failed." 00385 << "\nInput stream is probably mispositioned now." << std::endl; 00386 return is; 00387 } 00388 v.push_back(uu); 00389 } 00390 getState(v); 00391 return (is); 00392 } 00393 00394 // is >> theSeed; Removed, encompassed by possibleKeywordInput() 00395 char endMarker [MarkerLen]; 00396 for (int i=0; i<2; ++i) { 00397 is >> table[theSeed][i]; 00398 } 00399 is >> std::ws; 00400 is.width(MarkerLen); 00401 is >> endMarker; 00402 if (strcmp(endMarker,"RanecuEngine-end")) { 00403 is.clear(std::ios::badbit | is.rdstate()); 00404 std::cerr << "\nRanecuEngine state description incomplete." 00405 << "\nInput stream is probably mispositioned now." << std::endl; 00406 return is; 00407 } 00408 00409 seq = int(theSeed); 00410 return is; 00411 } 00412 00413 bool RanecuEngine::get (const std::vector<unsigned long> & v) { 00414 if ((v[0] & 0xffffffffUL) != engineIDulong<RanecuEngine>()) { 00415 std::cerr << 00416 "\nRanecuEngine get:state vector has wrong ID word - state unchanged\n"; 00417 return false; 00418 } 00419 return getState(v); 00420 } 00421 00422 bool RanecuEngine::getState (const std::vector<unsigned long> & v) { 00423 if (v.size() != VECTOR_STATE_SIZE ) { 00424 std::cerr << 00425 "\nRanecuEngine get:state vector has wrong length - state unchanged\n"; 00426 return false; 00427 } 00428 theSeed = v[1]; 00429 table[theSeed][0] = v[2]; 00430 table[theSeed][1] = v[3]; 00431 seq = int(theSeed); 00432 return true; 00433 } 00434 00435 00436 } // namespace CLHEP