CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: RanshiEngine.cc,v 1.4.4.2.2.1 2008/11/13 21:35:23 garren Exp $ 00002 // -*- C++ -*- 00003 // 00004 // ----------------------------------------------------------------------- 00005 // HEP Random 00006 // --- RanshiEngine --- 00007 // class implementation file 00008 // ----------------------------------------------------------------------- 00009 // 00010 // This algorithm implements the random number generator as proposed by 00011 // "F. Gutbrod, Comp. Phys. Comm. 87 (1995) 291-306". 00012 // 00013 // ======================================================================= 00014 // Ken Smith - Created: 9th June 1998 00015 // - Removed pow() from flat method: 21st Jul 1998 00016 // - Added conversion operators: 6th Aug 1998 00017 // J. Marraffino - Added some explicit casts to deal with 00018 // machines where sizeof(int) != sizeof(long) 22 Aug 1998 00019 // M. Fischler - Modified constructors taking seeds to not 00020 // depend on numEngines (same seeds should 00021 // produce same sequences). Default still 00022 // depends on numEngines. 16 Sep 1998 00023 // - Modified use of the various exponents of 2 00024 // to avoid per-instance space overhead and 00025 // correct the rounding procedure 16 Sep 1998 00026 // J. Marraffino - Remove dependence on hepString class 13 May 1999 00027 // M. Fischler - In restore, checkFile for file not found 03 Dec 2004 00028 // M. Fischler - Methods for instance save/restore 12/8/04 00029 // M. Fischler - split get() into tag validation and 00030 // getState() for anonymous restores 12/27/04 00031 // M. Fischler - State-saving using only ints, for portability 4/12/05 00032 // 00033 // ======================================================================= 00034 00035 #include "CLHEP/Random/defs.h" 00036 #include "CLHEP/Random/RanshiEngine.h" 00037 #include "CLHEP/Random/engineIDulong.h" 00038 #include <string.h> 00039 #include <cmath> // for ldexp() 00040 00041 using namespace std; 00042 00043 namespace CLHEP { 00044 00045 static const int MarkerLen = 64; // Enough room to hold a begin or end marker. 00046 00047 double RanshiEngine::twoToMinus_32; 00048 double RanshiEngine::twoToMinus_53; 00049 double RanshiEngine::nearlyTwoToMinus_54; 00050 00051 std::string RanshiEngine::name() const {return "RanshiEngine";} 00052 00053 void RanshiEngine::powersOfTwo() { 00054 twoToMinus_32 = ldexp (1.0, -32); 00055 twoToMinus_53 = ldexp (1.0, -53); 00056 nearlyTwoToMinus_54 = ldexp (1.0, -54) - ldexp (1.0, -100); 00057 } 00058 00059 // Number of instances with automatic seed selection 00060 int RanshiEngine::numEngines = 0; 00061 00062 RanshiEngine::RanshiEngine() : halfBuff(0), numFlats(0) 00063 { 00064 powersOfTwo(); 00065 int i = 0; 00066 while (i < numBuff) { 00067 buffer[i] = (unsigned int)(numEngines+19780503L*(i+1)); 00068 ++i; 00069 } 00070 theSeed = numEngines+19780503L*++i; 00071 redSpin = (unsigned int)(theSeed & 0xffffffff); 00072 ++numEngines; 00073 for( i = 0; i < 10000; ++i) flat(); // Warm-up by running thorugh 10000 nums 00074 } 00075 00076 RanshiEngine::RanshiEngine(std::istream& is) : halfBuff(0), numFlats(0) 00077 { 00078 is >> *this; 00079 } 00080 00081 RanshiEngine::RanshiEngine(long seed) : halfBuff(0), numFlats(0) 00082 { 00083 powersOfTwo(); 00084 for (int i = 0; i < numBuff; ++i) { 00085 buffer[i] = (unsigned int)seed&0xffffffff; 00086 } 00087 theSeed = seed; 00088 redSpin = (unsigned int)(theSeed & 0xffffffff); 00089 int j; 00090 for (j = 0; j < numBuff*20; ++j) { // "warm-up" for engine to hit 00091 flat(); // every ball on average 20X. 00092 } 00093 } 00094 00095 RanshiEngine::RanshiEngine(int rowIndex, int colIndex) 00096 : halfBuff(0), numFlats(0) 00097 { 00098 powersOfTwo(); 00099 int i = 0; 00100 while( i < numBuff ) { 00101 buffer[i] = (unsigned int)((rowIndex + (i+1)*(colIndex+8))&0xffffffff); 00102 ++i; 00103 } 00104 theSeed = rowIndex; 00105 redSpin = colIndex & 0xffffffff; 00106 for( i = 0; i < 100; ++i) flat(); // Warm-up by running thorugh 100 nums 00107 } 00108 00109 RanshiEngine::~RanshiEngine() { } 00110 00111 RanshiEngine::RanshiEngine( const RanshiEngine & p ) 00112 : halfBuff(0), numFlats(0) 00113 { 00114 *this = p; 00115 } 00116 00117 RanshiEngine & RanshiEngine::operator=( const RanshiEngine & p ) { 00118 if( this != &p ) { 00119 halfBuff = p.halfBuff; 00120 numFlats = p.numFlats; 00121 redSpin = p.redSpin; 00122 for( int i =0; i < numBuff; ++i) { 00123 buffer[i] = p.buffer[i]; 00124 } 00125 } 00126 return *this; 00127 } 00128 00129 double RanshiEngine::flat() { 00130 unsigned int redAngle = (((numBuff/2) - 1) & redSpin) + halfBuff; 00131 unsigned int blkSpin = buffer[redAngle] & 0xffffffff; 00132 unsigned int boostResult = blkSpin ^ redSpin; 00133 00134 buffer[redAngle] = ((blkSpin << 17) | (blkSpin >> (32-17))) ^ redSpin; 00135 00136 redSpin = (blkSpin + numFlats++) & 0xffffffff; 00137 halfBuff = numBuff/2 - halfBuff; 00138 00139 return ( blkSpin * twoToMinus_32 + // most significant part 00140 (boostResult>>11) * twoToMinus_53 + // fill in remaining bits 00141 nearlyTwoToMinus_54); // non-zero 00142 } 00143 00144 void RanshiEngine::flatArray(const int size, double* vect) { 00145 for (int i = 0; i < size; ++i) { 00146 vect[i] = flat(); 00147 } 00148 } 00149 00150 void RanshiEngine::setSeed(long seed, int) { 00151 *this = RanshiEngine(seed); 00152 } 00153 00154 void RanshiEngine::setSeeds(const long* seeds, int) { 00155 if (*seeds) { 00156 int i = 0; 00157 while (seeds[i] && i < numBuff) { 00158 buffer[i] = (unsigned int)seeds[i]; 00159 ++i; 00160 } 00161 while (i < numBuff) { 00162 buffer[i] = buffer[i-1]; 00163 ++i; 00164 } 00165 theSeed = seeds[0]; 00166 redSpin = (unsigned int)theSeed; 00167 } 00168 theSeeds = seeds; 00169 } 00170 00171 void RanshiEngine::saveStatus(const char filename[]) const { 00172 std::ofstream outFile(filename, std::ios::out); 00173 if (!outFile.bad()) { 00174 outFile << "Uvec\n"; 00175 std::vector<unsigned long> v = put(); 00176 #ifdef TRACE_IO 00177 std::cout << "Result of v = put() is:\n"; 00178 #endif 00179 for (unsigned int i=0; i<v.size(); ++i) { 00180 outFile << v[i] << "\n"; 00181 #ifdef TRACE_IO 00182 std::cout << v[i] << " "; 00183 if (i%6==0) std::cout << "\n"; 00184 #endif 00185 } 00186 #ifdef TRACE_IO 00187 std::cout << "\n"; 00188 #endif 00189 } 00190 #ifdef REMOVED 00191 if (!outFile.bad()) { 00192 outFile << std::setprecision(20) << theSeed << std::endl; 00193 for (int i = 0; i < numBuff; ++i) { 00194 outFile << buffer[i] << " "; 00195 } 00196 outFile << redSpin << " " << numFlats << " " << halfBuff << std::endl; 00197 } 00198 #endif 00199 } 00200 00201 void RanshiEngine::restoreStatus(const char filename[]) { 00202 std::ifstream inFile(filename, std::ios::in); 00203 if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) { 00204 std::cerr << " -- Engine state remains unchanged\n"; 00205 return; 00206 } 00207 if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) { 00208 std::vector<unsigned long> v; 00209 unsigned long xin; 00210 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) { 00211 inFile >> xin; 00212 #ifdef TRACE_IO 00213 std::cout << "ivec = " << ivec << " xin = " << xin << " "; 00214 if (ivec%3 == 0) std::cout << "\n"; 00215 #endif 00216 if (!inFile) { 00217 inFile.clear(std::ios::badbit | inFile.rdstate()); 00218 std::cerr << "\nRanshiEngine state (vector) description improper." 00219 << "\nrestoreStatus has failed." 00220 << "\nInput stream is probably mispositioned now." << std::endl; 00221 return; 00222 } 00223 v.push_back(xin); 00224 } 00225 getState(v); 00226 return; 00227 } 00228 00229 if (!inFile.bad()) { 00230 // inFile >> theSeed; removed -- encompased by possibleKeywordInput 00231 for (int i = 0; i < numBuff; ++i) { 00232 inFile >> buffer[i]; 00233 } 00234 inFile >> redSpin >> numFlats >> halfBuff; 00235 } 00236 } 00237 00238 void RanshiEngine::showStatus() const { 00239 std::cout << std::setprecision(20) << std::endl; 00240 std::cout << "----------- Ranshi engine status ----------" << std::endl; 00241 std::cout << "Initial seed = " << theSeed << std::endl; 00242 std::cout << "Current red spin = " << redSpin << std::endl; 00243 std::cout << "Values produced = " << numFlats << std::endl; 00244 std::cout << "Side of buffer = " << (halfBuff ? "upper" : "lower") 00245 << std::endl; 00246 std::cout << "Current buffer = " << std::endl; 00247 for (int i = 0; i < numBuff; i+=4) { 00248 std::cout << std::setw(10) << std::setiosflags(std::ios::right) 00249 << buffer[i] << std::setw(11) << buffer[i+1] << std::setw(11) 00250 << buffer[i+2] << std::setw(11) << buffer[i+3] << std::endl; 00251 } 00252 std::cout << "-------------------------------------------" << std::endl; 00253 } 00254 00255 RanshiEngine::operator float() { 00256 unsigned int redAngle = (((numBuff/2) - 1) & redSpin) + halfBuff; 00257 unsigned int blkSpin = buffer[redAngle] & 0xffffffff; 00258 00259 buffer[redAngle] = ((blkSpin << 17) | (blkSpin >> (32-17))) ^ redSpin; 00260 00261 redSpin = (blkSpin + numFlats++) & 0xffffffff; 00262 halfBuff = numBuff/2 - halfBuff; 00263 00264 return float(blkSpin * twoToMinus_32); 00265 } 00266 00267 RanshiEngine::operator unsigned int() { 00268 unsigned int redAngle = (((numBuff/2) - 1) & redSpin) + halfBuff; 00269 unsigned int blkSpin = buffer[redAngle] & 0xffffffff; 00270 00271 buffer[redAngle] = ((blkSpin << 17) | (blkSpin >> (32-17))) ^ redSpin; 00272 00273 redSpin = (blkSpin + numFlats++) & 0xffffffff; 00274 halfBuff = numBuff/2 - halfBuff; 00275 00276 return blkSpin; 00277 } 00278 00279 std::ostream& RanshiEngine::put (std::ostream& os ) const { 00280 char beginMarker[] = "RanshiEngine-begin"; 00281 os << beginMarker << "\nUvec\n"; 00282 std::vector<unsigned long> v = put(); 00283 for (unsigned int i=0; i<v.size(); ++i) { 00284 os << v[i] << "\n"; 00285 } 00286 return os; 00287 #ifdef REMOVED 00288 char endMarker[] = "RanshiEngine-end"; 00289 int pr=os.precision(20); 00290 os << " " << beginMarker << " "; 00291 00292 os << theSeed << "\n"; 00293 for (int i = 0; i < numBuff; ++i) { 00294 os << buffer[i] << "\n"; 00295 } 00296 os << redSpin << " " << numFlats << "\n" << halfBuff; 00297 00298 os << " " << endMarker << "\n"; 00299 os.precision(pr); 00300 return os; 00301 #endif 00302 } 00303 00304 std::vector<unsigned long> RanshiEngine::put () const { 00305 std::vector<unsigned long> v; 00306 v.push_back (engineIDulong<RanshiEngine>()); 00307 for (int i = 0; i < numBuff; ++i) { 00308 v.push_back(static_cast<unsigned long>(buffer[i])); 00309 } 00310 v.push_back(static_cast<unsigned long>(redSpin)); 00311 v.push_back(static_cast<unsigned long>(numFlats)); 00312 v.push_back(static_cast<unsigned long>(halfBuff)); 00313 return v; 00314 } 00315 00316 std::istream& RanshiEngine::get (std::istream& is) { 00317 char beginMarker [MarkerLen]; 00318 is >> std::ws; 00319 is.width(MarkerLen); // causes the next read to the char* to be <= 00320 // that many bytes, INCLUDING A TERMINATION \0 00321 // (Stroustrup, section 21.3.2) 00322 is >> beginMarker; 00323 if (strcmp(beginMarker,"RanshiEngine-begin")) { 00324 is.clear(std::ios::badbit | is.rdstate()); 00325 std::cerr << "\nInput mispositioned or" 00326 << "\nRanshiEngine state description missing or" 00327 << "\nwrong engine type found." << std::endl; 00328 return is; 00329 } 00330 return getState(is); 00331 } 00332 00333 std::string RanshiEngine::beginTag ( ) { 00334 return "RanshiEngine-begin"; 00335 } 00336 00337 std::istream& RanshiEngine::getState (std::istream& is) { 00338 if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) { 00339 std::vector<unsigned long> v; 00340 unsigned long uu; 00341 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) { 00342 is >> uu; 00343 if (!is) { 00344 is.clear(std::ios::badbit | is.rdstate()); 00345 std::cerr << "\nRanshiEngine state (vector) description improper." 00346 << "\ngetState() has failed." 00347 << "\nInput stream is probably mispositioned now." << std::endl; 00348 return is; 00349 } 00350 v.push_back(uu); 00351 } 00352 getState(v); 00353 return (is); 00354 } 00355 00356 // is >> theSeed; Removed, encompassed by possibleKeywordInput() 00357 00358 char endMarker [MarkerLen]; 00359 for (int i = 0; i < numBuff; ++i) { 00360 is >> buffer[i]; 00361 } 00362 is >> redSpin >> numFlats >> halfBuff; 00363 is >> std::ws; 00364 is.width(MarkerLen); 00365 is >> endMarker; 00366 if (strcmp(endMarker,"RanshiEngine-end")) { 00367 is.clear(std::ios::badbit | is.rdstate()); 00368 std::cerr << "\nRanshiEngine state description incomplete." 00369 << "\nInput stream is probably mispositioned now." << std::endl; 00370 return is; 00371 } 00372 return is; 00373 } 00374 00375 bool RanshiEngine::get (const std::vector<unsigned long> & v) { 00376 if ((v[0] & 0xffffffffUL) != engineIDulong<RanshiEngine>()) { 00377 std::cerr << 00378 "\nRanshiEngine get:state vector has wrong ID word - state unchanged\n"; 00379 return false; 00380 } 00381 return getState(v); 00382 } 00383 00384 bool RanshiEngine::getState (const std::vector<unsigned long> & v) { 00385 if (v.size() != VECTOR_STATE_SIZE ) { 00386 std::cerr << 00387 "\nRanshiEngine get:state vector has wrong length - state unchanged\n"; 00388 return false; 00389 } 00390 for (int i = 0; i < numBuff; ++i) { 00391 buffer[i] = v[i+1]; 00392 } 00393 redSpin = v[numBuff+1]; 00394 numFlats = v[numBuff+2]; 00395 halfBuff = v[numBuff+3]; 00396 return true; 00397 } 00398 00399 } // namespace CLHEP