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