CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: JamesRandom.cc,v 1.4.4.3.2.1 2008/11/13 21:35:23 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> 00045 #include <cmath> 00046 #include <stdlib.h> 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 { 00064 setSeed(seed,0); 00065 setSeeds(&theSeed,0); 00066 } 00067 00068 HepJamesRandom::HepJamesRandom() // 15 Feb. 1998 JMM 00069 { 00070 long seeds[2]; 00071 long seed; 00072 int cycle,curIndex; 00073 00074 cycle = abs(int(numEngines/maxIndex)); 00075 curIndex = abs(int(numEngines%maxIndex)); 00076 numEngines += 1; 00077 long mask = ((cycle & 0x007fffff) << 8); 00078 HepRandom::getTheTableSeeds( seeds, curIndex ); 00079 seed = seeds[0]^mask; 00080 setSeed(seed,0); 00081 setSeeds(&theSeed,0); 00082 } 00083 00084 HepJamesRandom::HepJamesRandom(int rowIndex, int colIndex) // 15 Feb. 1998 JMM 00085 { 00086 long seed; 00087 long seeds[2]; 00088 00089 int cycle = abs(int(rowIndex/maxIndex)); 00090 int row = abs(int(rowIndex%maxIndex)); 00091 int col = abs(int(colIndex%2)); 00092 long mask = ((cycle & 0x000007ff) << 20); 00093 HepRandom::getTheTableSeeds( seeds, row ); 00094 seed = (seeds[col])^mask; 00095 setSeed(seed,0); 00096 setSeeds(&theSeed,0); 00097 } 00098 00099 HepJamesRandom::HepJamesRandom(std::istream& is) 00100 { 00101 is >> *this; 00102 } 00103 00104 HepJamesRandom::~HepJamesRandom() {} 00105 00106 HepJamesRandom::HepJamesRandom(const HepJamesRandom &p) 00107 { 00108 int ipos, jpos; 00109 if ((this != &p) && (&p)) { 00110 theSeed = p.getSeed(); 00111 setSeeds(&theSeed,0); 00112 for (int i=0; i<97; ++i) 00113 u[i] = p.u[i]; 00114 c = p.c; cd = p.cd; cm = p.cm; 00115 jpos = p.j97; 00116 ipos = (64+jpos)%97; 00117 i97 = ipos; j97 = jpos; 00118 } 00119 } 00120 00121 HepJamesRandom & HepJamesRandom::operator = (const HepJamesRandom &p) 00122 { 00123 int ipos, jpos; 00124 if ((this != &p) && (&p)) { 00125 theSeed = p.getSeed(); 00126 setSeeds(&theSeed,0); 00127 for (int i=0; i<97; ++i) 00128 u[i] = p.u[i]; 00129 c = p.c; cd = p.cd; cm = p.cm; 00130 jpos = p.j97; 00131 ipos = (64+jpos)%97; 00132 i97 = ipos; j97 = jpos; 00133 } 00134 return *this; 00135 } 00136 00137 void HepJamesRandom::saveStatus( const char filename[] ) const 00138 { 00139 std::ofstream outFile( filename, std::ios::out ) ; 00140 00141 if (!outFile.bad()) { 00142 outFile << "Uvec\n"; 00143 std::vector<unsigned long> v = put(); 00144 #ifdef TRACE_IO 00145 std::cout << "Result of v = put() is:\n"; 00146 #endif 00147 for (unsigned int i=0; i<v.size(); ++i) { 00148 outFile << v[i] << "\n"; 00149 #ifdef TRACE_IO 00150 std::cout << v[i] << " "; 00151 if (i%6==0) std::cout << "\n"; 00152 #endif 00153 } 00154 #ifdef TRACE_IO 00155 std::cout << "\n"; 00156 #endif 00157 } 00158 #ifdef REMOVED 00159 int pos = j97; 00160 outFile << theSeed << std::endl; 00161 for (int i=0; i<97; ++i) 00162 outFile << std::setprecision(20) << u[i] << " "; 00163 outFile << std::endl; 00164 outFile << std::setprecision(20) << c << " "; 00165 outFile << std::setprecision(20) << cd << " "; 00166 outFile << std::setprecision(20) << cm << std::endl; 00167 outFile << pos << std::endl; 00168 #endif 00169 } 00170 00171 void HepJamesRandom::restoreStatus( const char filename[] ) 00172 { 00173 int ipos, jpos; 00174 std::ifstream inFile( filename, std::ios::in); 00175 if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) { 00176 std::cerr << " -- Engine state remains unchanged\n"; 00177 return; 00178 } 00179 if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) { 00180 std::vector<unsigned long> v; 00181 unsigned long xin; 00182 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) { 00183 inFile >> xin; 00184 #ifdef TRACE_IO 00185 std::cout << "ivec = " << ivec << " xin = " << xin << " "; 00186 if (ivec%3 == 0) std::cout << "\n"; 00187 #endif 00188 if (!inFile) { 00189 inFile.clear(std::ios::badbit | inFile.rdstate()); 00190 std::cerr << "\nJamesRandom state (vector) description improper." 00191 << "\nrestoreStatus has failed." 00192 << "\nInput stream is probably mispositioned now." << std::endl; 00193 return; 00194 } 00195 v.push_back(xin); 00196 } 00197 getState(v); 00198 return; 00199 } 00200 00201 if (!inFile.bad() && !inFile.eof()) { 00202 // inFile >> theSeed; removed -- encompased by possibleKeywordInput 00203 for (int i=0; i<97; ++i) 00204 inFile >> u[i]; 00205 inFile >> c; inFile >> cd; inFile >> cm; 00206 inFile >> jpos; 00207 ipos = (64+jpos)%97; 00208 i97 = ipos; 00209 j97 = jpos; 00210 } 00211 } 00212 00213 void HepJamesRandom::showStatus() const 00214 { 00215 std::cout << std::endl; 00216 std::cout << "----- HepJamesRandom engine status -----" << std::endl; 00217 std::cout << " Initial seed = " << theSeed << std::endl; 00218 std::cout << " u[] = "; 00219 for (int i=0; i<97; ++i) 00220 std::cout << u[i] << " "; 00221 std::cout << std::endl; 00222 std::cout << " c = " << c << ", cd = " << cd << ", cm = " << cm 00223 << std::endl; 00224 std::cout << " i97 = " << i97 << ", u[i97] = " << u[i97] << std::endl; 00225 std::cout << " j97 = " << j97 << ", u[j97] = " << u[j97] << std::endl; 00226 std::cout << "----------------------------------------" << std::endl; 00227 } 00228 00229 void HepJamesRandom::setSeed(long seed, int) 00230 { 00231 // The input value for "seed" should be within the range [0,900000000] 00232 // 00233 // Negative seeds result in serious flaws in the randomness; 00234 // seeds above 900000000 are OK because of the %177 in the expression for i, 00235 // but may have the same effect as other seeds below 900000000. 00236 00237 int m, n; 00238 float s, t; 00239 long mm; 00240 00241 if (seed < 0) { 00242 std::cout << "Seed for HepJamesRandom must be non-negative\n" 00243 << "Seed value supplied was " << seed 00244 << "\nUsing its absolute value instead\n"; 00245 seed = -seed; 00246 } 00247 00248 long ij = seed/30082; 00249 long kl = seed - 30082*ij; 00250 long i = (ij/177) % 177 + 2; 00251 long j = ij % 177 + 2; 00252 long k = (kl/169) % 178 + 1; 00253 long l = kl % 169; 00254 00255 theSeed = seed; 00256 00257 for ( n = 1 ; n < 98 ; n++ ) { 00258 s = 0.0; 00259 t = 0.5; 00260 for ( m = 1 ; m < 25 ; m++) { 00261 mm = ( ( (i*j) % 179 ) * k ) % 179; 00262 i = j; 00263 j = k; 00264 k = mm; 00265 l = ( 53 * l + 1 ) % 169; 00266 if ( (l*mm % 64 ) >= 32 ) 00267 s += t; 00268 t *= 0.5; 00269 } 00270 u[n-1] = s; 00271 } 00272 c = 362436.0 / 16777216.0; 00273 cd = 7654321.0 / 16777216.0; 00274 cm = 16777213.0 / 16777216.0; 00275 00276 i97 = 96; 00277 j97 = 32; 00278 00279 } 00280 00281 void HepJamesRandom::setSeeds(const long* seeds, int) 00282 { 00283 setSeed(seeds ? *seeds : 19780503L, 0); 00284 theSeeds = seeds; 00285 } 00286 00287 double HepJamesRandom::flat() 00288 { 00289 double uni; 00290 00291 do { 00292 uni = u[i97] - u[j97]; 00293 if ( uni < 0.0 ) uni++; 00294 u[i97] = uni; 00295 00296 if (i97 == 0) i97 = 96; 00297 else i97--; 00298 00299 if (j97 == 0) j97 = 96; 00300 else j97--; 00301 00302 c -= cd; 00303 if (c < 0.0) c += cm; 00304 00305 uni -= c; 00306 if (uni < 0.0) uni += 1.0; 00307 } while ( uni <= 0.0 || uni >= 1.0 ); 00308 00309 return uni; 00310 } 00311 00312 void HepJamesRandom::flatArray(const int size, double* vect) 00313 { 00314 // double uni; 00315 int i; 00316 00317 for (i=0; i<size; ++i) { 00318 vect[i] = flat(); 00319 } 00320 } 00321 00322 HepJamesRandom::operator unsigned int() { 00323 return ((unsigned int)(flat() * exponent_bit_32) & 0xffffffff ) | 00324 (((unsigned int)( u[i97] * exponent_bit_32)>>16) & 0xff); 00325 } 00326 00327 std::ostream & HepJamesRandom::put ( std::ostream& os ) const { 00328 char beginMarker[] = "JamesRandom-begin"; 00329 os << beginMarker << "\nUvec\n"; 00330 std::vector<unsigned long> v = put(); 00331 for (unsigned int i=0; i<v.size(); ++i) { 00332 os << v[i] << "\n"; 00333 } 00334 return os; 00335 #ifdef REMOVED 00336 char endMarker[] = "JamesRandom-end"; 00337 int pos = j97; 00338 int pr = os.precision(20); 00339 os << " " << beginMarker << " "; 00340 os << theSeed << " "; 00341 for (int i=0; i<97; ++i) { 00342 os << std::setprecision(20) << u[i] << "\n"; 00343 } 00344 os << std::setprecision(20) << c << " "; 00345 os << std::setprecision(20) << cd << " "; 00346 os << std::setprecision(20) << cm << " "; 00347 os << pos << "\n"; 00348 os << endMarker << "\n"; 00349 os.precision(pr); 00350 return os; 00351 #endif 00352 } 00353 00354 std::vector<unsigned long> HepJamesRandom::put () const { 00355 std::vector<unsigned long> v; 00356 v.push_back (engineIDulong<HepJamesRandom>()); 00357 std::vector<unsigned long> t; 00358 for (int i=0; i<97; ++i) { 00359 t = DoubConv::dto2longs(u[i]); 00360 v.push_back(t[0]); v.push_back(t[1]); 00361 } 00362 t = DoubConv::dto2longs(c); 00363 v.push_back(t[0]); v.push_back(t[1]); 00364 t = DoubConv::dto2longs(cd); 00365 v.push_back(t[0]); v.push_back(t[1]); 00366 t = DoubConv::dto2longs(cm); 00367 v.push_back(t[0]); v.push_back(t[1]); 00368 v.push_back(static_cast<unsigned long>(j97)); 00369 return v; 00370 } 00371 00372 00373 std::istream & HepJamesRandom::get ( std::istream& is) { 00374 char beginMarker [MarkerLen]; 00375 is >> std::ws; 00376 is.width(MarkerLen); // causes the next read to the char* to be <= 00377 // that many bytes, INCLUDING A TERMINATION \0 00378 // (Stroustrup, section 21.3.2) 00379 is >> beginMarker; 00380 if (strcmp(beginMarker,"JamesRandom-begin")) { 00381 is.clear(std::ios::badbit | is.rdstate()); 00382 std::cerr << "\nInput stream mispositioned or" 00383 << "\nJamesRandom state description missing or" 00384 << "\nwrong engine type found." << std::endl; 00385 return is; 00386 } 00387 return getState(is); 00388 } 00389 00390 std::string HepJamesRandom::beginTag ( ) { 00391 return "JamesRandom-begin"; 00392 } 00393 00394 std::istream & HepJamesRandom::getState ( std::istream& is) { 00395 if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) { 00396 std::vector<unsigned long> v; 00397 unsigned long uu; 00398 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) { 00399 is >> uu; 00400 if (!is) { 00401 is.clear(std::ios::badbit | is.rdstate()); 00402 std::cerr << "\nJamesRandom state (vector) description improper." 00403 << "\ngetState() has failed." 00404 << "\nInput stream is probably mispositioned now." << std::endl; 00405 return is; 00406 } 00407 v.push_back(uu); 00408 } 00409 getState(v); 00410 return (is); 00411 } 00412 00413 // is >> theSeed; Removed, encompassed by possibleKeywordInput() 00414 00415 int ipos, jpos; 00416 char endMarker [MarkerLen]; 00417 for (int i=0; i<97; ++i) { 00418 is >> u[i]; 00419 } 00420 is >> c; is >> cd; is >> cm; 00421 is >> jpos; 00422 is >> std::ws; 00423 is.width(MarkerLen); 00424 is >> endMarker; 00425 if(strcmp(endMarker,"JamesRandom-end")) { 00426 is.clear(std::ios::badbit | is.rdstate()); 00427 std::cerr << "\nJamesRandom state description incomplete." 00428 << "\nInput stream is probably mispositioned now." << std::endl; 00429 return is; 00430 } 00431 00432 ipos = (64+jpos)%97; 00433 i97 = ipos; 00434 j97 = jpos; 00435 return is; 00436 } 00437 00438 bool HepJamesRandom::get (const std::vector<unsigned long> & v) { 00439 if ( (v[0] & 0xffffffffUL) != engineIDulong<HepJamesRandom>()) { 00440 std::cerr << 00441 "\nHepJamesRandom get:state vector has wrong ID word - state unchanged\n"; 00442 return false; 00443 } 00444 return getState(v); 00445 } 00446 00447 bool HepJamesRandom::getState (const std::vector<unsigned long> & v) { 00448 if (v.size() != VECTOR_STATE_SIZE ) { 00449 std::cerr << 00450 "\nHepJamesRandom get:state vector has wrong length - state unchanged\n"; 00451 return false; 00452 } 00453 std::vector<unsigned long> t(2); 00454 for (int i=0; i<97; ++i) { 00455 t[0] = v[2*i+1]; t[1] = v[2*i+2]; 00456 u[i] = DoubConv::longs2double(t); 00457 } 00458 t[0] = v[195]; t[1] = v[196]; c = DoubConv::longs2double(t); 00459 t[0] = v[197]; t[1] = v[198]; cd = DoubConv::longs2double(t); 00460 t[0] = v[199]; t[1] = v[200]; cm = DoubConv::longs2double(t); 00461 j97 = v[201]; 00462 i97 = (64+j97)%97; 00463 return true; 00464 } 00465 00466 } // namespace CLHEP