CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: JamesRandom.cc,v 1.6 2010/06/16 17:24:53 garren Exp $ 00002 // -*- C++ -*- 00003 // 00004 // ----------------------------------------------------------------------- 00005 // HEP Random 00006 // --- HepJamesRandom --- 00007 // class implementation file 00008 // ----------------------------------------------------------------------- 00009 // This file is part of Geant4 (simulation toolkit for HEP). 00010 // 00011 // This algorithm implements the original universal random number generator 00012 // as proposed by Marsaglia & Zaman in report FSU-SCRI-87-50 and coded 00013 // in FORTRAN77 by Fred James as the RANMAR generator, part of the MATHLIB 00014 // HEP library. 00015 00016 // ======================================================================= 00017 // Gabriele Cosmo - Created: 5th September 1995 00018 // - Fixed a bug in setSeed(): 26th February 1996 00019 // - Minor corrections: 31st October 1996 00020 // - Added methods for engine status: 19th November 1996 00021 // - Fixed bug in setSeeds(): 15th September 1997 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 // V. Innocente - changed pointers to indices 3 may 2000 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 - Enforcement that seeds be non-negative 00033 // (lest the sequence be non-random) 2/14/05 00034 // M. Fischler - put/get for vectors of ulongs 3/14/05 00035 // M. Fischler - State-saving using only ints, for portability 4/12/05 00036 // 00037 // ======================================================================= 00038 00039 #include "CLHEP/Random/defs.h" 00040 #include "CLHEP/Random/Random.h" 00041 #include "CLHEP/Random/JamesRandom.h" 00042 #include "CLHEP/Random/engineIDulong.h" 00043 #include "CLHEP/Random/DoubConv.hh" 00044 #include <string.h> // for strcmp 00045 #include <cmath> 00046 #include <cstdlib> 00047 00048 //#define TRACE_IO 00049 00050 namespace CLHEP { 00051 00052 static const int MarkerLen = 64; // Enough room to hold a begin or end marker. 00053 00054 std::string HepJamesRandom::name() const {return "HepJamesRandom";} 00055 00056 // Number of instances with automatic seed selection 00057 int HepJamesRandom::numEngines = 0; 00058 00059 // Maximum index into the seed table 00060 int HepJamesRandom::maxIndex = 215; 00061 00062 HepJamesRandom::HepJamesRandom(long seed) 00063 : HepRandomEngine() 00064 { 00065 setSeed(seed,0); 00066 setSeeds(&theSeed,0); 00067 } 00068 00069 HepJamesRandom::HepJamesRandom() // 15 Feb. 1998 JMM 00070 : HepRandomEngine() 00071 { 00072 long seeds[2]; 00073 long seed; 00074 00075 int cycle = std::abs(int(numEngines/maxIndex)); 00076 int curIndex = std::abs(int(numEngines%maxIndex)); 00077 ++numEngines; 00078 long mask = ((cycle & 0x007fffff) << 8); 00079 HepRandom::getTheTableSeeds( seeds, curIndex ); 00080 seed = seeds[0]^mask; 00081 setSeed(seed,0); 00082 setSeeds(&theSeed,0); 00083 } 00084 00085 HepJamesRandom::HepJamesRandom(int rowIndex, int colIndex) // 15 Feb. 1998 JMM 00086 : HepRandomEngine() 00087 { 00088 long seed; 00089 long seeds[2]; 00090 00091 int cycle = std::abs(int(rowIndex/maxIndex)); 00092 int row = std::abs(int(rowIndex%maxIndex)); 00093 int col = std::abs(int(colIndex%2)); 00094 long mask = ((cycle & 0x000007ff) << 20); 00095 HepRandom::getTheTableSeeds( seeds, row ); 00096 seed = (seeds[col])^mask; 00097 setSeed(seed,0); 00098 setSeeds(&theSeed,0); 00099 } 00100 00101 HepJamesRandom::HepJamesRandom(std::istream& is) 00102 : HepRandomEngine() 00103 { 00104 is >> *this; 00105 } 00106 00107 HepJamesRandom::~HepJamesRandom() {} 00108 00109 void HepJamesRandom::saveStatus( const char filename[] ) const 00110 { 00111 std::ofstream outFile( filename, std::ios::out ) ; 00112 00113 if (!outFile.bad()) { 00114 outFile << "Uvec\n"; 00115 std::vector<unsigned long> v = put(); 00116 #ifdef TRACE_IO 00117 std::cout << "Result of v = put() is:\n"; 00118 #endif 00119 for (unsigned int i=0; i<v.size(); ++i) { 00120 outFile << v[i] << "\n"; 00121 #ifdef TRACE_IO 00122 std::cout << v[i] << " "; 00123 if (i%6==0) std::cout << "\n"; 00124 #endif 00125 } 00126 #ifdef TRACE_IO 00127 std::cout << "\n"; 00128 #endif 00129 } 00130 #ifdef REMOVED 00131 int pos = j97; 00132 outFile << theSeed << std::endl; 00133 for (int i=0; i<97; ++i) 00134 outFile << std::setprecision(20) << u[i] << " "; 00135 outFile << std::endl; 00136 outFile << std::setprecision(20) << c << " "; 00137 outFile << std::setprecision(20) << cd << " "; 00138 outFile << std::setprecision(20) << cm << std::endl; 00139 outFile << pos << std::endl; 00140 #endif 00141 } 00142 00143 void HepJamesRandom::restoreStatus( const char filename[] ) 00144 { 00145 int ipos, jpos; 00146 std::ifstream inFile( filename, std::ios::in); 00147 if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) { 00148 std::cerr << " -- Engine state remains unchanged\n"; 00149 return; 00150 } 00151 if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) { 00152 std::vector<unsigned long> v; 00153 unsigned long xin; 00154 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) { 00155 inFile >> xin; 00156 #ifdef TRACE_IO 00157 std::cout << "ivec = " << ivec << " xin = " << xin << " "; 00158 if (ivec%3 == 0) std::cout << "\n"; 00159 #endif 00160 if (!inFile) { 00161 inFile.clear(std::ios::badbit | inFile.rdstate()); 00162 std::cerr << "\nJamesRandom state (vector) description improper." 00163 << "\nrestoreStatus has failed." 00164 << "\nInput stream is probably mispositioned now." << std::endl; 00165 return; 00166 } 00167 v.push_back(xin); 00168 } 00169 getState(v); 00170 return; 00171 } 00172 00173 if (!inFile.bad() && !inFile.eof()) { 00174 // inFile >> theSeed; removed -- encompased by possibleKeywordInput 00175 for (int i=0; i<97; ++i) 00176 inFile >> u[i]; 00177 inFile >> c; inFile >> cd; inFile >> cm; 00178 inFile >> jpos; 00179 ipos = (64+jpos)%97; 00180 i97 = ipos; 00181 j97 = jpos; 00182 } 00183 } 00184 00185 void HepJamesRandom::showStatus() const 00186 { 00187 std::cout << std::endl; 00188 std::cout << "----- HepJamesRandom engine status -----" << std::endl; 00189 std::cout << " Initial seed = " << theSeed << std::endl; 00190 std::cout << " u[] = "; 00191 for (int i=0; i<97; ++i) 00192 std::cout << u[i] << " "; 00193 std::cout << std::endl; 00194 std::cout << " c = " << c << ", cd = " << cd << ", cm = " << cm 00195 << std::endl; 00196 std::cout << " i97 = " << i97 << ", u[i97] = " << u[i97] << std::endl; 00197 std::cout << " j97 = " << j97 << ", u[j97] = " << u[j97] << std::endl; 00198 std::cout << "----------------------------------------" << std::endl; 00199 } 00200 00201 void HepJamesRandom::setSeed(long seed, int) 00202 { 00203 // The input value for "seed" should be within the range [0,900000000] 00204 // 00205 // Negative seeds result in serious flaws in the randomness; 00206 // seeds above 900000000 are OK because of the %177 in the expression for i, 00207 // but may have the same effect as other seeds below 900000000. 00208 00209 int m, n; 00210 float s, t; 00211 long mm; 00212 00213 if (seed < 0) { 00214 std::cout << "Seed for HepJamesRandom must be non-negative\n" 00215 << "Seed value supplied was " << seed 00216 << "\nUsing its absolute value instead\n"; 00217 seed = -seed; 00218 } 00219 00220 long ij = seed/30082; 00221 long kl = seed - 30082*ij; 00222 long i = (ij/177) % 177 + 2; 00223 long j = ij % 177 + 2; 00224 long k = (kl/169) % 178 + 1; 00225 long l = kl % 169; 00226 00227 theSeed = seed; 00228 00229 for ( n = 1 ; n < 98 ; n++ ) { 00230 s = 0.0; 00231 t = 0.5; 00232 for ( m = 1 ; m < 25 ; m++) { 00233 mm = ( ( (i*j) % 179 ) * k ) % 179; 00234 i = j; 00235 j = k; 00236 k = mm; 00237 l = ( 53 * l + 1 ) % 169; 00238 if ( (l*mm % 64 ) >= 32 ) 00239 s += t; 00240 t *= 0.5; 00241 } 00242 u[n-1] = s; 00243 } 00244 c = 362436.0 / 16777216.0; 00245 cd = 7654321.0 / 16777216.0; 00246 cm = 16777213.0 / 16777216.0; 00247 00248 i97 = 96; 00249 j97 = 32; 00250 00251 } 00252 00253 void HepJamesRandom::setSeeds(const long* seeds, int) 00254 { 00255 setSeed(seeds ? *seeds : 19780503L, 0); 00256 theSeeds = seeds; 00257 } 00258 00259 double HepJamesRandom::flat() 00260 { 00261 double uni; 00262 00263 do { 00264 uni = u[i97] - u[j97]; 00265 if ( uni < 0.0 ) uni++; 00266 u[i97] = uni; 00267 00268 if (i97 == 0) i97 = 96; 00269 else i97--; 00270 00271 if (j97 == 0) j97 = 96; 00272 else j97--; 00273 00274 c -= cd; 00275 if (c < 0.0) c += cm; 00276 00277 uni -= c; 00278 if (uni < 0.0) uni += 1.0; 00279 } while ( uni <= 0.0 || uni >= 1.0 ); 00280 00281 return uni; 00282 } 00283 00284 void HepJamesRandom::flatArray(const int size, double* vect) 00285 { 00286 // double uni; 00287 int i; 00288 00289 for (i=0; i<size; ++i) { 00290 vect[i] = flat(); 00291 } 00292 } 00293 00294 HepJamesRandom::operator unsigned int() { 00295 return ((unsigned int)(flat() * exponent_bit_32()) & 0xffffffff ) | 00296 (((unsigned int)( u[i97] * exponent_bit_32())>>16) & 0xff); 00297 } 00298 00299 std::ostream & HepJamesRandom::put ( std::ostream& os ) const { 00300 char beginMarker[] = "JamesRandom-begin"; 00301 os << beginMarker << "\nUvec\n"; 00302 std::vector<unsigned long> v = put(); 00303 for (unsigned int i=0; i<v.size(); ++i) { 00304 os << v[i] << "\n"; 00305 } 00306 return os; 00307 #ifdef REMOVED 00308 char endMarker[] = "JamesRandom-end"; 00309 int pos = j97; 00310 int pr = os.precision(20); 00311 os << " " << beginMarker << " "; 00312 os << theSeed << " "; 00313 for (int i=0; i<97; ++i) { 00314 os << std::setprecision(20) << u[i] << "\n"; 00315 } 00316 os << std::setprecision(20) << c << " "; 00317 os << std::setprecision(20) << cd << " "; 00318 os << std::setprecision(20) << cm << " "; 00319 os << pos << "\n"; 00320 os << endMarker << "\n"; 00321 os.precision(pr); 00322 return os; 00323 #endif 00324 } 00325 00326 std::vector<unsigned long> HepJamesRandom::put () const { 00327 std::vector<unsigned long> v; 00328 v.push_back (engineIDulong<HepJamesRandom>()); 00329 std::vector<unsigned long> t; 00330 for (int i=0; i<97; ++i) { 00331 t = DoubConv::dto2longs(u[i]); 00332 v.push_back(t[0]); v.push_back(t[1]); 00333 } 00334 t = DoubConv::dto2longs(c); 00335 v.push_back(t[0]); v.push_back(t[1]); 00336 t = DoubConv::dto2longs(cd); 00337 v.push_back(t[0]); v.push_back(t[1]); 00338 t = DoubConv::dto2longs(cm); 00339 v.push_back(t[0]); v.push_back(t[1]); 00340 v.push_back(static_cast<unsigned long>(j97)); 00341 return v; 00342 } 00343 00344 00345 std::istream & HepJamesRandom::get ( std::istream& is) { 00346 char beginMarker [MarkerLen]; 00347 is >> std::ws; 00348 is.width(MarkerLen); // causes the next read to the char* to be <= 00349 // that many bytes, INCLUDING A TERMINATION \0 00350 // (Stroustrup, section 21.3.2) 00351 is >> beginMarker; 00352 if (strcmp(beginMarker,"JamesRandom-begin")) { 00353 is.clear(std::ios::badbit | is.rdstate()); 00354 std::cerr << "\nInput stream mispositioned or" 00355 << "\nJamesRandom state description missing or" 00356 << "\nwrong engine type found." << std::endl; 00357 return is; 00358 } 00359 return getState(is); 00360 } 00361 00362 std::string HepJamesRandom::beginTag ( ) { 00363 return "JamesRandom-begin"; 00364 } 00365 00366 std::istream & HepJamesRandom::getState ( std::istream& is) { 00367 if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) { 00368 std::vector<unsigned long> v; 00369 unsigned long uu; 00370 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) { 00371 is >> uu; 00372 if (!is) { 00373 is.clear(std::ios::badbit | is.rdstate()); 00374 std::cerr << "\nJamesRandom state (vector) description improper." 00375 << "\ngetState() has failed." 00376 << "\nInput stream is probably mispositioned now." << std::endl; 00377 return is; 00378 } 00379 v.push_back(uu); 00380 } 00381 getState(v); 00382 return (is); 00383 } 00384 00385 // is >> theSeed; Removed, encompassed by possibleKeywordInput() 00386 00387 int ipos, jpos; 00388 char endMarker [MarkerLen]; 00389 for (int i=0; i<97; ++i) { 00390 is >> u[i]; 00391 } 00392 is >> c; is >> cd; is >> cm; 00393 is >> jpos; 00394 is >> std::ws; 00395 is.width(MarkerLen); 00396 is >> endMarker; 00397 if(strcmp(endMarker,"JamesRandom-end")) { 00398 is.clear(std::ios::badbit | is.rdstate()); 00399 std::cerr << "\nJamesRandom state description incomplete." 00400 << "\nInput stream is probably mispositioned now." << std::endl; 00401 return is; 00402 } 00403 00404 ipos = (64+jpos)%97; 00405 i97 = ipos; 00406 j97 = jpos; 00407 return is; 00408 } 00409 00410 bool HepJamesRandom::get (const std::vector<unsigned long> & v) { 00411 if ( (v[0] & 0xffffffffUL) != engineIDulong<HepJamesRandom>()) { 00412 std::cerr << 00413 "\nHepJamesRandom get:state vector has wrong ID word - state unchanged\n"; 00414 return false; 00415 } 00416 return getState(v); 00417 } 00418 00419 bool HepJamesRandom::getState (const std::vector<unsigned long> & v) { 00420 if (v.size() != VECTOR_STATE_SIZE ) { 00421 std::cerr << 00422 "\nHepJamesRandom get:state vector has wrong length - state unchanged\n"; 00423 return false; 00424 } 00425 std::vector<unsigned long> t(2); 00426 for (int i=0; i<97; ++i) { 00427 t[0] = v[2*i+1]; t[1] = v[2*i+2]; 00428 u[i] = DoubConv::longs2double(t); 00429 } 00430 t[0] = v[195]; t[1] = v[196]; c = DoubConv::longs2double(t); 00431 t[0] = v[197]; t[1] = v[198]; cd = DoubConv::longs2double(t); 00432 t[0] = v[199]; t[1] = v[200]; cm = DoubConv::longs2double(t); 00433 j97 = v[201]; 00434 i97 = (64+j97)%97; 00435 return true; 00436 } 00437 00438 } // namespace CLHEP