CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: Ranlux64Engine.cc,v 1.7 2010/10/21 21:32:02 garren Exp $ 00002 // -*- C++ -*- 00003 // 00004 // ----------------------------------------------------------------------- 00005 // HEP Random 00006 // --- Ranlux64Engine --- 00007 // class implementation file 00008 // ----------------------------------------------------------------------- 00009 // A double-precision implementation of the RanluxEngine generator as 00010 // decsribed by the notes of the original ranlux author (Martin Luscher) 00011 // 00012 // See the note by Martin Luscher, December 1997, entitiled 00013 // Double-precision implementation of the random number generator ranlux 00014 // 00015 // ======================================================================= 00016 // Ken Smith - Initial draft: 14th Jul 1998 00017 // - Removed std::pow() from flat method 14th Jul 1998 00018 // - Added conversion operators: 6th Aug 1998 00019 // 00020 // Mark Fischler The following were modified mostly to make the routine 00021 // exactly match the Luscher algorithm in generating 48-bit 00022 // randoms: 00023 // 9/9/98 - Substantial changes in what used to be flat() to match 00024 // algorithm in Luscher's ranlxd.c 00025 // - Added update() method for 12 numbers, making flat() trivial 00026 // - Added advance() method to hold the unrolled loop for update 00027 // - Distinction between three forms of seeding such that it 00028 // is impossible to get same sequence from different forms - 00029 // done by discarding some fraction of one macro cycle which 00030 // is different for the three cases 00031 // - Change the misnomer "seed_table" to the more accurate 00032 // "randoms" 00033 // - Removed the no longer needed count12, i_lag, j_lag, etc. 00034 // - Corrected seed procedure which had been filling bits past 00035 // 2^-48. This actually was very bad, invalidating the 00036 // number theory behind the proof that ranlxd is good. 00037 // - Addition of 2**(-49) to generated number to prevent zero 00038 // from being returned; this does not affect the sequence 00039 // itself. 00040 // - Corrected ecu seeding, which had been supplying only 00041 // numbers less than 1/2. This is probably moot. 00042 // 9/15/98 - Modified use of the various exponents of 2 00043 // to avoid per-instance space overhead. Note that these 00044 // are initialized in setSeed, which EVERY constructor 00045 // must invoke. 00046 // J. Marraffino - Remove dependence on hepString class 13 May 1999 00047 // M. Fischler - In restore, checkFile for file not found 03 Dec 2004 00048 // M. Fischler - put get Methods for distrib instance save/restore 12/8/04 00049 // M. Fischler - split get() into tag validation and 00050 // getState() for anonymous restores 12/27/04 00051 // M. Fischler - put/get for vectors of ulongs 3/14/05 00052 // M. Fischler - State-saving using only ints, for portability 4/12/05 00053 // 00054 // ======================================================================= 00055 00056 #include "CLHEP/Random/defs.h" 00057 #include "CLHEP/Random/Random.h" 00058 #include "CLHEP/Random/Ranlux64Engine.h" 00059 #include "CLHEP/Random/engineIDulong.h" 00060 #include "CLHEP/Random/DoubConv.hh" 00061 #include <string.h> // for strcmp 00062 #include <cstdlib> // for std::abs(int) 00063 #include <limits> // for numeric_limits 00064 00065 namespace CLHEP { 00066 00067 static const int MarkerLen = 64; // Enough room to hold a begin or end marker. 00068 00069 00070 // Number of instances with automatic seed selection 00071 int Ranlux64Engine::numEngines = 0; 00072 00073 // Maximum index into the seed table 00074 int Ranlux64Engine::maxIndex = 215; 00075 00076 #ifndef USING_VISUAL 00077 namespace detail { 00078 00079 template< std::size_t n, 00080 bool = n < std::size_t(std::numeric_limits<unsigned long>::digits) > 00081 struct do_right_shift; 00082 template< std::size_t n > 00083 struct do_right_shift<n,true> 00084 { 00085 unsigned long operator()(unsigned long value) { return value >> n; } 00086 }; 00087 template< std::size_t n > 00088 struct do_right_shift<n,false> 00089 { 00090 unsigned long operator()(unsigned long) { return 0ul; } 00091 }; 00092 00093 template< std::size_t nbits > 00094 unsigned long rshift( unsigned long value ) 00095 { return do_right_shift<nbits>()(value); } 00096 00097 } // namespace detail 00098 #endif 00099 00100 std::string Ranlux64Engine::name() const {return "Ranlux64Engine";} 00101 00102 Ranlux64Engine::Ranlux64Engine() 00103 : HepRandomEngine() 00104 { 00105 luxury = 1; 00106 int cycle = std::abs(int(numEngines/maxIndex)); 00107 int curIndex = std::abs(int(numEngines%maxIndex)); 00108 numEngines +=1; 00109 long mask = ((cycle & 0x007fffff) << 8); 00110 long seedlist[2]; 00111 HepRandom::getTheTableSeeds( seedlist, curIndex ); 00112 seedlist[0] ^= mask; 00113 seedlist[1] = 0; 00114 00115 setSeeds(seedlist, luxury); 00116 advance ( 8 ); // Discard some iterations and ensure that 00117 // this sequence won't match one where seeds 00118 // were provided. 00119 } 00120 00121 Ranlux64Engine::Ranlux64Engine(long seed, int lux) 00122 : HepRandomEngine() 00123 { 00124 luxury = lux; 00125 long seedlist[2]={seed,0}; 00126 setSeeds(seedlist, lux); 00127 advance ( 2*lux + 1 ); // Discard some iterations to use a different 00128 // point in the sequence. 00129 } 00130 00131 Ranlux64Engine::Ranlux64Engine(int rowIndex, int, int lux) 00132 : HepRandomEngine() 00133 { 00134 luxury = lux; 00135 int cycle = std::abs(int(rowIndex/maxIndex)); 00136 int row = std::abs(int(rowIndex%maxIndex)); 00137 long mask = (( cycle & 0x000007ff ) << 20 ); 00138 long seedlist[2]; 00139 HepRandom::getTheTableSeeds( seedlist, row ); 00140 seedlist[0] ^= mask; 00141 seedlist[1]= 0; 00142 setSeeds(seedlist, lux); 00143 } 00144 00145 Ranlux64Engine::Ranlux64Engine( std::istream& is ) 00146 : HepRandomEngine() 00147 { 00148 is >> *this; 00149 } 00150 00151 Ranlux64Engine::~Ranlux64Engine() {} 00152 00153 double Ranlux64Engine::flat() { 00154 // Luscher improves the speed by computing several numbers in a shot, 00155 // in a manner similar to that of the Tausworth in DualRand or the Hurd 00156 // engines. Thus, the real work is done in update(). Here we merely ensure 00157 // that zero, which the algorithm can produce, is never returned by flat(). 00158 00159 if (index <= 0) update(); 00160 return randoms[--index] + twoToMinus_49(); 00161 } 00162 00163 void Ranlux64Engine::update() { 00164 // Update the stash of twelve random numbers. 00165 // When this routione is entered, index is always 0. The randoms 00166 // contains the last 12 numbers in the sequents: s[0] is x[a+11], 00167 // s[1] is x[a+10] ... and s[11] is x[a] for some a. Carry contains 00168 // the last carry value (c[a+11]). 00169 // 00170 // The recursion relation (3) in Luscher's note says 00171 // delta[n] = x[n-s] = x[n-r] -c[n-1] or for n=a+12, 00172 // delta[a+12] = x[a+7] - x[a] -c[a+11] where we use r=12, s=5 per eqn. (7) 00173 // This reduces to 00174 // s[11] = s[4] - s[11] - carry. 00175 // The next number similarly will be given by s[10] = s[3] - s[10] - carry, 00176 // and so forth until s[0] is filled. 00177 // 00178 // However, we need to skip 397, 202 or 109 numbers - these are not divisible 00179 // by 12 - to "fare well in the spectral test". 00180 00181 advance(pDozens); 00182 00183 // Since we wish at the end to have the 12 last numbers in the order of 00184 // s[11] first, till s[0] last, we will have to do 1, 10, or 1 iterations 00185 // and then re-arrange to place to get the oldest one in s[11]. 00186 // Generically, this will imply re-arranging the s array at the end, 00187 // but we can treat the special case of endIters = 1 separately for superior 00188 // efficiency in the cases of levels 0 and 2. 00189 00190 register double y1; 00191 00192 if ( endIters == 1 ) { // Luxury levels 0 and 2 will go here 00193 y1 = randoms[ 4] - randoms[11] - carry; 00194 if ( y1 < 0.0 ) { 00195 y1 += 1.0; 00196 carry = twoToMinus_48(); 00197 } else { 00198 carry = 0.0; 00199 } 00200 randoms[11] = randoms[10]; 00201 randoms[10] = randoms[ 9]; 00202 randoms[ 9] = randoms[ 8]; 00203 randoms[ 8] = randoms[ 7]; 00204 randoms[ 7] = randoms[ 6]; 00205 randoms[ 6] = randoms[ 5]; 00206 randoms[ 5] = randoms[ 4]; 00207 randoms[ 4] = randoms[ 3]; 00208 randoms[ 3] = randoms[ 2]; 00209 randoms[ 2] = randoms[ 1]; 00210 randoms[ 1] = randoms[ 0]; 00211 randoms[ 0] = y1; 00212 00213 } else { 00214 00215 int m, nr, ns; 00216 for ( m = 0, nr = 11, ns = 4; m < endIters; ++m, --nr ) { 00217 y1 = randoms [ns] - randoms[nr] - carry; 00218 if ( y1 < 0.0 ) { 00219 y1 += 1.0; 00220 carry = twoToMinus_48(); 00221 } else { 00222 carry = 0.0; 00223 } 00224 randoms[nr] = y1; 00225 --ns; 00226 if ( ns < 0 ) { 00227 ns = 11; 00228 } 00229 } // loop on m 00230 00231 double temp[12]; 00232 for (m=0; m<12; m++) { 00233 temp[m]=randoms[m]; 00234 } 00235 00236 ns = 11 - endIters; 00237 for (m=11; m>=0; --m) { 00238 randoms[m] = temp[ns]; 00239 --ns; 00240 if ( ns < 0 ) { 00241 ns = 11; 00242 } 00243 } 00244 00245 } 00246 00247 // Now when we return, there are 12 fresh usable numbers in s[11] ... s[0] 00248 00249 index = 11; 00250 00251 } // update() 00252 00253 void Ranlux64Engine::advance(int dozens) { 00254 00255 register double y1, y2, y3; 00256 register double cValue = twoToMinus_48(); 00257 register double zero = 0.0; 00258 register double one = 1.0; 00259 00260 // Technical note: We use Luscher's trick to only do the 00261 // carry subtraction when we really have to. Like him, we use 00262 // three registers instead of two so that we avoid sequences 00263 // like storing y1 then immediately replacing its value: 00264 // some architectures lose time when this is done. 00265 00266 // Luscher's ranlxd.c fills the stash going 00267 // upward. We fill it downward to save a bit of time in the 00268 // flat() routine at no cost later. This means that while 00269 // Luscher's ir is jr+5, our n-r is (n-s)-5. (Note that 00270 // though ranlxd.c initializes ir and jr to 11 and 7, ir as 00271 // used is 5 more than jr because update is entered after 00272 // incrementing ir.) 00273 // 00274 00275 // I have CAREFULLY checked that the algorithms do match 00276 // in all details. 00277 00278 int k; 00279 for ( k = dozens; k > 0; --k ) { 00280 00281 y1 = randoms[ 4] - randoms[11] - carry; 00282 00283 y2 = randoms[ 3] - randoms[10]; 00284 if ( y1 < zero ) { 00285 y1 += one; 00286 y2 -= cValue; 00287 } 00288 randoms[11] = y1; 00289 00290 y3 = randoms[ 2] - randoms[ 9]; 00291 if ( y2 < zero ) { 00292 y2 += one; 00293 y3 -= cValue; 00294 } 00295 randoms[10] = y2; 00296 00297 y1 = randoms[ 1] - randoms[ 8]; 00298 if ( y3 < zero ) { 00299 y3 += one; 00300 y1 -= cValue; 00301 } 00302 randoms[ 9] = y3; 00303 00304 y2 = randoms[ 0] - randoms[ 7]; 00305 if ( y1 < zero ) { 00306 y1 += one; 00307 y2 -= cValue; 00308 } 00309 randoms[ 8] = y1; 00310 00311 y3 = randoms[11] - randoms[ 6]; 00312 if ( y2 < zero ) { 00313 y2 += one; 00314 y3 -= cValue; 00315 } 00316 randoms[ 7] = y2; 00317 00318 y1 = randoms[10] - randoms[ 5]; 00319 if ( y3 < zero ) { 00320 y3 += one; 00321 y1 -= cValue; 00322 } 00323 randoms[ 6] = y3; 00324 00325 y2 = randoms[ 9] - randoms[ 4]; 00326 if ( y1 < zero ) { 00327 y1 += one; 00328 y2 -= cValue; 00329 } 00330 randoms[ 5] = y1; 00331 00332 y3 = randoms[ 8] - randoms[ 3]; 00333 if ( y2 < zero ) { 00334 y2 += one; 00335 y3 -= cValue; 00336 } 00337 randoms[ 4] = y2; 00338 00339 y1 = randoms[ 7] - randoms[ 2]; 00340 if ( y3 < zero ) { 00341 y3 += one; 00342 y1 -= cValue; 00343 } 00344 randoms[ 3] = y3; 00345 00346 y2 = randoms[ 6] - randoms[ 1]; 00347 if ( y1 < zero ) { 00348 y1 += one; 00349 y2 -= cValue; 00350 } 00351 randoms[ 2] = y1; 00352 00353 y3 = randoms[ 5] - randoms[ 0]; 00354 if ( y2 < zero ) { 00355 y2 += one; 00356 y3 -= cValue; 00357 } 00358 randoms[ 1] = y2; 00359 00360 if ( y3 < zero ) { 00361 y3 += one; 00362 carry = cValue; 00363 } 00364 randoms[ 0] = y3; 00365 00366 } // End of major k loop doing 12 numbers at each cycle 00367 00368 } // advance(dozens) 00369 00370 void Ranlux64Engine::flatArray(const int size, double* vect) { 00371 for( int i=0; i < size; ++i ) { 00372 vect[i] = flat(); 00373 } 00374 } 00375 00376 void Ranlux64Engine::setSeed(long seed, int lux) { 00377 00378 // The initialization is carried out using a Multiplicative 00379 // Congruential generator using formula constants of L'Ecuyer 00380 // as described in "A review of pseudorandom number generators" 00381 // (Fred James) published in Computer Physics Communications 60 (1990) 00382 // pages 329-344 00383 00384 const int ecuyer_a(53668); 00385 const int ecuyer_b(40014); 00386 const int ecuyer_c(12211); 00387 const int ecuyer_d(2147483563); 00388 00389 const int lux_levels[3] = {109, 202, 397}; 00390 theSeed = seed; 00391 00392 if( (lux > 2)||(lux < 0) ){ 00393 pDiscard = (lux >= 12) ? (lux-12) : lux_levels[1]; 00394 }else{ 00395 pDiscard = lux_levels[luxury]; 00396 } 00397 pDozens = pDiscard / 12; 00398 endIters = pDiscard % 12; 00399 00400 long init_table[24]; 00401 long next_seed = seed; 00402 long k_multiple; 00403 int i; 00404 next_seed &= 0xffffffff; 00405 while( next_seed >= ecuyer_d ) { 00406 next_seed -= ecuyer_d; 00407 } 00408 00409 for(i = 0;i != 24;i++){ 00410 k_multiple = next_seed / ecuyer_a; 00411 next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a) 00412 - k_multiple * ecuyer_c; 00413 if(next_seed < 0) { 00414 next_seed += ecuyer_d; 00415 } 00416 next_seed &= 0xffffffff; 00417 init_table[i] = next_seed; 00418 } 00419 // are we on a 64bit machine? 00420 if( sizeof(long) >= 8 ) { 00421 long topbits1, topbits2; 00422 #ifdef USING_VISUAL 00423 topbits1 = ( seed >> 32) & 0xffff ; 00424 topbits2 = ( seed >> 48) & 0xffff ; 00425 #else 00426 topbits1 = detail::rshift<32>(seed) & 0xffff ; 00427 topbits2 = detail::rshift<48>(seed) & 0xffff ; 00428 #endif 00429 init_table[0] ^= topbits1; 00430 init_table[2] ^= topbits2; 00431 //std::cout << " init_table[0] " << init_table[0] << " from " << topbits1 << std::endl; 00432 //std::cout << " init_table[2] " << init_table[2] << " from " << topbits2 << std::endl; 00433 } 00434 00435 for(i = 0;i < 12; i++){ 00436 randoms[i] = (init_table[2*i ] ) * 2.0 * twoToMinus_32() + 00437 (init_table[2*i+1] >> 15) * twoToMinus_48(); 00438 //if( randoms[i] < 0. || randoms[i] > 1. ) { 00439 //std::cout << "setSeed: init_table " << init_table[2*i ] << std::endl; 00440 //std::cout << "setSeed: init_table " << init_table[2*i+1] << std::endl; 00441 //std::cout << "setSeed: random " << i << " is " << randoms[i] << std::endl; 00442 //} 00443 } 00444 00445 carry = 0.0; 00446 if ( randoms[11] == 0. ) carry = twoToMinus_48(); 00447 index = 11; 00448 00449 } // setSeed() 00450 00451 void Ranlux64Engine::setSeeds(const long * seeds, int lux) { 00452 // old code only uses the first long in seeds 00453 // setSeed( *seeds ? *seeds : 32767, lux ); 00454 // theSeeds = seeds; 00455 00456 // using code from Ranlux - even those are 32bit seeds, 00457 // that is good enough to completely differentiate the sequences 00458 00459 const int ecuyer_a = 53668; 00460 const int ecuyer_b = 40014; 00461 const int ecuyer_c = 12211; 00462 const int ecuyer_d = 2147483563; 00463 00464 const int lux_levels[3] = {109, 202, 397}; 00465 const long *seedptr; 00466 00467 theSeeds = seeds; 00468 seedptr = seeds; 00469 00470 if(seeds == 0){ 00471 setSeed(theSeed,lux); 00472 theSeeds = &theSeed; 00473 return; 00474 } 00475 00476 theSeed = *seeds; 00477 00478 // number of additional random numbers that need to be 'thrown away' 00479 // every 24 numbers is set using luxury level variable. 00480 00481 if( (lux > 2)||(lux < 0) ){ 00482 pDiscard = (lux >= 12) ? (lux-12) : lux_levels[1]; 00483 }else{ 00484 pDiscard = lux_levels[luxury]; 00485 } 00486 pDozens = pDiscard / 12; 00487 endIters = pDiscard % 12; 00488 00489 long init_table[24]; 00490 long next_seed = *seeds; 00491 long k_multiple; 00492 int i; 00493 00494 for( i = 0;(i != 24)&&(*seedptr != 0);i++){ 00495 init_table[i] = *seedptr & 0xffffffff; 00496 seedptr++; 00497 } 00498 00499 if(i != 24){ 00500 next_seed = init_table[i-1]; 00501 for(;i != 24;i++){ 00502 k_multiple = next_seed / ecuyer_a; 00503 next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a) 00504 - k_multiple * ecuyer_c; 00505 if(next_seed < 0) { 00506 next_seed += ecuyer_d; 00507 } 00508 next_seed &= 0xffffffff; 00509 init_table[i] = next_seed; 00510 } 00511 } 00512 00513 for(i = 0;i < 12; i++){ 00514 randoms[i] = (init_table[2*i ] ) * 2.0 * twoToMinus_32() + 00515 (init_table[2*i+1] >> 15) * twoToMinus_48(); 00516 } 00517 00518 carry = 0.0; 00519 if ( randoms[11] == 0. ) carry = twoToMinus_48(); 00520 index = 11; 00521 00522 } 00523 00524 void Ranlux64Engine::saveStatus( const char filename[] ) const 00525 { 00526 std::ofstream outFile( filename, std::ios::out ) ; 00527 if (!outFile.bad()) { 00528 outFile << "Uvec\n"; 00529 std::vector<unsigned long> v = put(); 00530 #ifdef TRACE_IO 00531 std::cout << "Result of v = put() is:\n"; 00532 #endif 00533 for (unsigned int i=0; i<v.size(); ++i) { 00534 outFile << v[i] << "\n"; 00535 #ifdef TRACE_IO 00536 std::cout << v[i] << " "; 00537 if (i%6==0) std::cout << "\n"; 00538 #endif 00539 } 00540 #ifdef TRACE_IO 00541 std::cout << "\n"; 00542 #endif 00543 } 00544 #ifdef REMOVED 00545 if (!outFile.bad()) { 00546 outFile << theSeed << std::endl; 00547 for (int i=0; i<12; ++i) { 00548 outFile <<std::setprecision(20) << randoms[i] << std::endl; 00549 } 00550 outFile << std::setprecision(20) << carry << " " << index << std::endl; 00551 outFile << luxury << " " << pDiscard << std::endl; 00552 } 00553 #endif 00554 } 00555 00556 void Ranlux64Engine::restoreStatus( const char filename[] ) 00557 { 00558 std::ifstream inFile( filename, std::ios::in); 00559 if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) { 00560 std::cerr << " -- Engine state remains unchanged\n"; 00561 return; 00562 } 00563 if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) { 00564 std::vector<unsigned long> v; 00565 unsigned long xin; 00566 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) { 00567 inFile >> xin; 00568 #ifdef TRACE_IO 00569 std::cout << "ivec = " << ivec << " xin = " << xin << " "; 00570 if (ivec%3 == 0) std::cout << "\n"; 00571 #endif 00572 if (!inFile) { 00573 inFile.clear(std::ios::badbit | inFile.rdstate()); 00574 std::cerr << "\nJamesRandom state (vector) description improper." 00575 << "\nrestoreStatus has failed." 00576 << "\nInput stream is probably mispositioned now." << std::endl; 00577 return; 00578 } 00579 v.push_back(xin); 00580 } 00581 getState(v); 00582 return; 00583 } 00584 00585 if (!inFile.bad() && !inFile.eof()) { 00586 // inFile >> theSeed; removed -- encompased by possibleKeywordInput 00587 for (int i=0; i<12; ++i) { 00588 inFile >> randoms[i]; 00589 } 00590 inFile >> carry; inFile >> index; 00591 inFile >> luxury; inFile >> pDiscard; 00592 pDozens = pDiscard / 12; 00593 endIters = pDiscard % 12; 00594 } 00595 } 00596 00597 void Ranlux64Engine::showStatus() const 00598 { 00599 std::cout << std::endl; 00600 std::cout << "--------- Ranlux engine status ---------" << std::endl; 00601 std::cout << " Initial seed = " << theSeed << std::endl; 00602 std::cout << " randoms[] = "; 00603 for (int i=0; i<12; ++i) { 00604 std::cout << randoms[i] << std::endl; 00605 } 00606 std::cout << std::endl; 00607 std::cout << " carry = " << carry << ", index = " << index << std::endl; 00608 std::cout << " luxury = " << luxury << " pDiscard = " 00609 << pDiscard << std::endl; 00610 std::cout << "----------------------------------------" << std::endl; 00611 } 00612 00613 std::ostream & Ranlux64Engine::put( std::ostream& os ) const 00614 { 00615 char beginMarker[] = "Ranlux64Engine-begin"; 00616 os << beginMarker << "\nUvec\n"; 00617 std::vector<unsigned long> v = put(); 00618 for (unsigned int i=0; i<v.size(); ++i) { 00619 os << v[i] << "\n"; 00620 } 00621 return os; 00622 #ifdef REMOVED 00623 char endMarker[] = "Ranlux64Engine-end"; 00624 int pr = os.precision(20); 00625 os << " " << beginMarker << " "; 00626 os << theSeed << " "; 00627 for (int i=0; i<12; ++i) { 00628 os << randoms[i] << std::endl; 00629 } 00630 os << carry << " " << index << " "; 00631 os << luxury << " " << pDiscard << "\n"; 00632 os << endMarker << " "; 00633 os.precision(pr); 00634 return os; 00635 #endif 00636 } 00637 00638 std::vector<unsigned long> Ranlux64Engine::put () const { 00639 std::vector<unsigned long> v; 00640 v.push_back (engineIDulong<Ranlux64Engine>()); 00641 std::vector<unsigned long> t; 00642 for (int i=0; i<12; ++i) { 00643 t = DoubConv::dto2longs(randoms[i]); 00644 v.push_back(t[0]); v.push_back(t[1]); 00645 } 00646 t = DoubConv::dto2longs(carry); 00647 v.push_back(t[0]); v.push_back(t[1]); 00648 v.push_back(static_cast<unsigned long>(index)); 00649 v.push_back(static_cast<unsigned long>(luxury)); 00650 v.push_back(static_cast<unsigned long>(pDiscard)); 00651 return v; 00652 } 00653 00654 std::istream & Ranlux64Engine::get ( std::istream& is ) 00655 { 00656 char beginMarker [MarkerLen]; 00657 is >> std::ws; 00658 is.width(MarkerLen); // causes the next read to the char* to be <= 00659 // that many bytes, INCLUDING A TERMINATION \0 00660 // (Stroustrup, section 21.3.2) 00661 is >> beginMarker; 00662 if (strcmp(beginMarker,"Ranlux64Engine-begin")) { 00663 is.clear(std::ios::badbit | is.rdstate()); 00664 std::cerr << "\nInput stream mispositioned or" 00665 << "\nRanlux64Engine state description missing or" 00666 << "\nwrong engine type found." << std::endl; 00667 return is; 00668 } 00669 return getState(is); 00670 } 00671 00672 std::string Ranlux64Engine::beginTag ( ) { 00673 return "Ranlux64Engine-begin"; 00674 } 00675 00676 std::istream & Ranlux64Engine::getState ( std::istream& is ) 00677 { 00678 if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) { 00679 std::vector<unsigned long> v; 00680 unsigned long uu; 00681 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) { 00682 is >> uu; 00683 if (!is) { 00684 is.clear(std::ios::badbit | is.rdstate()); 00685 std::cerr << "\nRanlux64Engine state (vector) description improper." 00686 << "\ngetState() has failed." 00687 << "\nInput stream is probably mispositioned now." << std::endl; 00688 return is; 00689 } 00690 v.push_back(uu); 00691 } 00692 getState(v); 00693 return (is); 00694 } 00695 00696 // is >> theSeed; Removed, encompassed by possibleKeywordInput() 00697 00698 char endMarker [MarkerLen]; 00699 for (int i=0; i<12; ++i) { 00700 is >> randoms[i]; 00701 } 00702 is >> carry; is >> index; 00703 is >> luxury; is >> pDiscard; 00704 pDozens = pDiscard / 12; 00705 endIters = pDiscard % 12; 00706 is >> std::ws; 00707 is.width(MarkerLen); 00708 is >> endMarker; 00709 if (strcmp(endMarker,"Ranlux64Engine-end")) { 00710 is.clear(std::ios::badbit | is.rdstate()); 00711 std::cerr << "\nRanlux64Engine state description incomplete." 00712 << "\nInput stream is probably mispositioned now." << std::endl; 00713 return is; 00714 } 00715 return is; 00716 } 00717 00718 bool Ranlux64Engine::get (const std::vector<unsigned long> & v) { 00719 if ((v[0] & 0xffffffffUL) != engineIDulong<Ranlux64Engine>()) { 00720 std::cerr << 00721 "\nRanlux64Engine get:state vector has wrong ID word - state unchanged\n"; 00722 return false; 00723 } 00724 return getState(v); 00725 } 00726 00727 bool Ranlux64Engine::getState (const std::vector<unsigned long> & v) { 00728 if (v.size() != VECTOR_STATE_SIZE ) { 00729 std::cerr << 00730 "\nRanlux64Engine get:state vector has wrong length - state unchanged\n"; 00731 return false; 00732 } 00733 std::vector<unsigned long> t(2); 00734 for (int i=0; i<12; ++i) { 00735 t[0] = v[2*i+1]; t[1] = v[2*i+2]; 00736 randoms[i] = DoubConv::longs2double(t); 00737 } 00738 t[0] = v[25]; t[1] = v[26]; 00739 carry = DoubConv::longs2double(t); 00740 index = v[27]; 00741 luxury = v[28]; 00742 pDiscard = v[29]; 00743 return true; 00744 } 00745 00746 } // namespace CLHEP