CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: RanecuEngine.cc,v 1.7 2010/07/20 18:06:02 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> // for strcmp 00045 #include <cmath> 00046 #include <cstdlib> 00047 00048 namespace CLHEP { 00049 00050 static const int MarkerLen = 64; // Enough room to hold a begin or end marker. 00051 00052 static const double prec = 4.6566128E-10; 00053 00054 std::string RanecuEngine::name() const {return "RanecuEngine";} 00055 00056 void RanecuEngine::further_randomize (int seq1, int col, int index, int modulus) 00057 { 00058 table[seq1][col] -= (index&0x3FFFFFFF); 00059 while (table[seq1][col] <= 0) table[seq1][col] += (modulus-1); 00060 } // mf 6/22/10 00061 00062 // Number of instances with automatic seed selection 00063 int RanecuEngine::numEngines = 0; 00064 00065 RanecuEngine::RanecuEngine() 00066 : HepRandomEngine() 00067 { 00068 int cycle = std::abs(int(numEngines/maxSeq)); 00069 seq = std::abs(int(numEngines%maxSeq)); 00070 numEngines += 1; 00071 theSeed = seq; 00072 long mask = ((cycle & 0x007fffff) << 8); 00073 for (int i=0; i<2; ++i) { 00074 for (int j=0; j<maxSeq; ++j) { 00075 HepRandom::getTheTableSeeds(table[j],j); 00076 table[j][i] ^= mask; 00077 } 00078 } 00079 theSeeds = &table[seq][0]; 00080 } 00081 00082 RanecuEngine::RanecuEngine(int index) 00083 : HepRandomEngine() 00084 { 00085 int cycle = std::abs(int(index/maxSeq)); 00086 seq = std::abs(int(index%maxSeq)); 00087 theSeed = seq; 00088 long mask = ((cycle & 0x000007ff) << 20); 00089 for (int j=0; j<maxSeq; ++j) { 00090 HepRandom::getTheTableSeeds(table[j],j); 00091 table[j][0] ^= mask; 00092 table[j][1] ^= mask; 00093 } 00094 theSeeds = &table[seq][0]; 00095 further_randomize (seq, 0, index, shift1); // mf 6/22/10 00096 } 00097 00098 RanecuEngine::RanecuEngine(std::istream& is) 00099 : HepRandomEngine() 00100 { 00101 is >> *this; 00102 } 00103 00104 RanecuEngine::~RanecuEngine() {} 00105 00106 void RanecuEngine::setSeed(long index, int dum) 00107 { 00108 seq = std::abs(int(index%maxSeq)); 00109 theSeed = seq; 00110 HepRandom::getTheTableSeeds(table[seq],seq); 00111 theSeeds = &table[seq][0]; 00112 further_randomize (seq, 0, index, shift1); // mf 6/22/10 00113 further_randomize (seq, 1, dum, shift2); // mf 6/22/10 00114 } 00115 00116 void RanecuEngine::setSeeds(const long* seeds, int pos) 00117 { 00118 if (pos != -1) { 00119 seq = std::abs(int(pos%maxSeq)); 00120 theSeed = seq; 00121 } 00122 // only positive seeds are allowed 00123 table[seq][0] = std::abs(seeds[0])%shift1; 00124 table[seq][1] = std::abs(seeds[1])%shift2; 00125 theSeeds = &table[seq][0]; 00126 } 00127 00128 void RanecuEngine::setIndex(long index) 00129 { 00130 seq = std::abs(int(index%maxSeq)); 00131 theSeed = seq; 00132 theSeeds = &table[seq][0]; 00133 } 00134 00135 void RanecuEngine::saveStatus( const char filename[] ) const 00136 { 00137 std::ofstream outFile( filename, std::ios::out ) ; 00138 00139 if (!outFile.bad()) { 00140 outFile << "Uvec\n"; 00141 std::vector<unsigned long> v = put(); 00142 #ifdef TRACE_IO 00143 std::cout << "Result of v = put() is:\n"; 00144 #endif 00145 for (unsigned int i=0; i<v.size(); ++i) { 00146 outFile << v[i] << "\n"; 00147 #ifdef TRACE_IO 00148 std::cout << v[i] << " "; 00149 if (i%6==0) std::cout << "\n"; 00150 #endif 00151 } 00152 #ifdef TRACE_IO 00153 std::cout << "\n"; 00154 #endif 00155 } 00156 #ifdef REMOVED 00157 if (!outFile.bad()) { 00158 outFile << theSeed << std::endl; 00159 for (int i=0; i<2; ++i) 00160 outFile << table[theSeed][i] << " "; 00161 outFile << std::endl; 00162 } 00163 #endif 00164 } 00165 00166 void RanecuEngine::restoreStatus( const char filename[] ) 00167 { 00168 std::ifstream inFile( filename, std::ios::in); 00169 if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) { 00170 std::cerr << " -- Engine state remains unchanged\n"; 00171 return; 00172 } 00173 if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) { 00174 std::vector<unsigned long> v; 00175 unsigned long xin; 00176 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) { 00177 inFile >> xin; 00178 #ifdef TRACE_IO 00179 std::cout << "ivec = " << ivec << " xin = " << xin << " "; 00180 if (ivec%3 == 0) std::cout << "\n"; 00181 #endif 00182 if (!inFile) { 00183 inFile.clear(std::ios::badbit | inFile.rdstate()); 00184 std::cerr << "\nJamesRandom state (vector) description improper." 00185 << "\nrestoreStatus has failed." 00186 << "\nInput stream is probably mispositioned now." << std::endl; 00187 return; 00188 } 00189 v.push_back(xin); 00190 } 00191 getState(v); 00192 return; 00193 } 00194 00195 if (!inFile.bad() && !inFile.eof()) { 00196 // inFile >> theSeed; removed -- encompased by possibleKeywordInput 00197 for (int i=0; i<2; ++i) 00198 inFile >> table[theSeed][i]; 00199 seq = int(theSeed); 00200 } 00201 } 00202 00203 void RanecuEngine::showStatus() const 00204 { 00205 std::cout << std::endl; 00206 std::cout << "--------- Ranecu engine status ---------" << std::endl; 00207 std::cout << " Initial seed (index) = " << theSeed << std::endl; 00208 std::cout << " Current couple of seeds = " 00209 << table[theSeed][0] << ", " 00210 << table[theSeed][1] << std::endl; 00211 std::cout << "----------------------------------------" << std::endl; 00212 } 00213 00214 double RanecuEngine::flat() 00215 { 00216 const int index = seq; 00217 long seed1 = table[index][0]; 00218 long seed2 = table[index][1]; 00219 00220 int k1 = (int)(seed1/ecuyer_b); 00221 int k2 = (int)(seed2/ecuyer_e); 00222 00223 seed1 = ecuyer_a*(seed1-k1*ecuyer_b)-k1*ecuyer_c; 00224 if (seed1 < 0) seed1 += shift1; 00225 seed2 = ecuyer_d*(seed2-k2*ecuyer_e)-k2*ecuyer_f; 00226 if (seed2 < 0) seed2 += shift2; 00227 00228 table[index][0] = seed1; 00229 table[index][1] = seed2; 00230 00231 long diff = seed1-seed2; 00232 00233 if (diff <= 0) diff += (shift1-1); 00234 return (double)(diff*prec); 00235 } 00236 00237 void RanecuEngine::flatArray(const int size, double* vect) 00238 { 00239 const int index = seq; 00240 long seed1 = table[index][0]; 00241 long seed2 = table[index][1]; 00242 int k1, k2; 00243 register int i; 00244 00245 for (i=0; i<size; ++i) 00246 { 00247 k1 = (int)(seed1/ecuyer_b); 00248 k2 = (int)(seed2/ecuyer_e); 00249 00250 seed1 = ecuyer_a*(seed1-k1*ecuyer_b)-k1*ecuyer_c; 00251 if (seed1 < 0) seed1 += shift1; 00252 seed2 = ecuyer_d*(seed2-k2*ecuyer_e)-k2*ecuyer_f; 00253 if (seed2 < 0) seed2 += shift2; 00254 00255 long diff = seed1-seed2; 00256 if (diff <= 0) diff += (shift1-1); 00257 00258 vect[i] = (double)(diff*prec); 00259 } 00260 table[index][0] = seed1; 00261 table[index][1] = seed2; 00262 } 00263 00264 RanecuEngine::operator unsigned int() { 00265 const int index = seq; 00266 long seed1 = table[index][0]; 00267 long seed2 = table[index][1]; 00268 00269 int k1 = (int)(seed1/ecuyer_b); 00270 int k2 = (int)(seed2/ecuyer_e); 00271 00272 seed1 = ecuyer_a*(seed1-k1*ecuyer_b)-k1*ecuyer_c; 00273 if (seed1 < 0) seed1 += shift1; 00274 seed2 = ecuyer_d*(seed2-k2*ecuyer_e)-k2*ecuyer_f; 00275 if (seed2 < 0) seed2 += shift2; 00276 00277 table[index][0] = seed1; 00278 table[index][1] = seed2; 00279 long diff = seed1-seed2; 00280 if( diff <= 0 ) diff += (shift1-1); 00281 00282 return ((diff << 1) | (seed1&1))& 0xffffffff; 00283 } 00284 00285 std::ostream & RanecuEngine::put( std::ostream& os ) const 00286 { 00287 char beginMarker[] = "RanecuEngine-begin"; 00288 os << beginMarker << "\nUvec\n"; 00289 std::vector<unsigned long> v = put(); 00290 for (unsigned int i=0; i<v.size(); ++i) { 00291 os << v[i] << "\n"; 00292 } 00293 return os; 00294 #ifdef REMOVED 00295 char endMarker[] = "RanecuEngine-end"; 00296 os << " " << beginMarker << "\n"; 00297 os << theSeed << " "; 00298 for (int i=0; i<2; ++i) { 00299 os << table[theSeed][i] << "\n"; 00300 } 00301 os << endMarker << "\n"; 00302 return os; 00303 #endif 00304 } 00305 00306 std::vector<unsigned long> RanecuEngine::put () const { 00307 std::vector<unsigned long> v; 00308 v.push_back (engineIDulong<RanecuEngine>()); 00309 v.push_back(static_cast<unsigned long>(theSeed)); 00310 v.push_back(static_cast<unsigned long>(table[theSeed][0])); 00311 v.push_back(static_cast<unsigned long>(table[theSeed][1])); 00312 return v; 00313 } 00314 00315 std::istream & RanecuEngine::get ( std::istream& is ) 00316 { 00317 char beginMarker [MarkerLen]; 00318 00319 is >> std::ws; 00320 is.width(MarkerLen); // causes the next read to the char* to be <= 00321 // that many bytes, INCLUDING A TERMINATION \0 00322 // (Stroustrup, section 21.3.2) 00323 is >> beginMarker; 00324 if (strcmp(beginMarker,"RanecuEngine-begin")) { 00325 is.clear(std::ios::badbit | is.rdstate()); 00326 std::cerr << "\nInput stream mispositioned or" 00327 << "\nRanecuEngine state description missing or" 00328 << "\nwrong engine type found." << std::endl; 00329 return is; 00330 } 00331 return getState(is); 00332 } 00333 00334 std::string RanecuEngine::beginTag ( ) { 00335 return "RanecuEngine-begin"; 00336 } 00337 00338 std::istream & RanecuEngine::getState ( std::istream& is ) 00339 { 00340 if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) { 00341 std::vector<unsigned long> v; 00342 unsigned long uu; 00343 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) { 00344 is >> uu; 00345 if (!is) { 00346 is.clear(std::ios::badbit | is.rdstate()); 00347 std::cerr << "\nRanecuEngine state (vector) description improper." 00348 << "\ngetState() has failed." 00349 << "\nInput stream is probably mispositioned now." << std::endl; 00350 return is; 00351 } 00352 v.push_back(uu); 00353 } 00354 getState(v); 00355 return (is); 00356 } 00357 00358 // is >> theSeed; Removed, encompassed by possibleKeywordInput() 00359 char endMarker [MarkerLen]; 00360 for (int i=0; i<2; ++i) { 00361 is >> table[theSeed][i]; 00362 } 00363 is >> std::ws; 00364 is.width(MarkerLen); 00365 is >> endMarker; 00366 if (strcmp(endMarker,"RanecuEngine-end")) { 00367 is.clear(std::ios::badbit | is.rdstate()); 00368 std::cerr << "\nRanecuEngine state description incomplete." 00369 << "\nInput stream is probably mispositioned now." << std::endl; 00370 return is; 00371 } 00372 00373 seq = int(theSeed); 00374 return is; 00375 } 00376 00377 bool RanecuEngine::get (const std::vector<unsigned long> & v) { 00378 if ((v[0] & 0xffffffffUL) != engineIDulong<RanecuEngine>()) { 00379 std::cerr << 00380 "\nRanecuEngine get:state vector has wrong ID word - state unchanged\n"; 00381 return false; 00382 } 00383 return getState(v); 00384 } 00385 00386 bool RanecuEngine::getState (const std::vector<unsigned long> & v) { 00387 if (v.size() != VECTOR_STATE_SIZE ) { 00388 std::cerr << 00389 "\nRanecuEngine get:state vector has wrong length - state unchanged\n"; 00390 return false; 00391 } 00392 theSeed = v[1]; 00393 table[theSeed][0] = v[2]; 00394 table[theSeed][1] = v[3]; 00395 seq = int(theSeed); 00396 return true; 00397 } 00398 00399 00400 } // namespace CLHEP