CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

Ranlux64Engine.cc
Go to the documentation of this file.
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