Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/externals/clhep/src/Ranlux64Engine.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /externals/clhep/src/Ranlux64Engine.cc (Version 11.3.0) and /externals/clhep/src/Ranlux64Engine.cc (Version 9.6)


                                                   >>   1 // $Id:$
  1 // -*- C++ -*-                                      2 // -*- C++ -*-
  2 //                                                  3 //
  3 // -------------------------------------------      4 // -----------------------------------------------------------------------
  4 //                             HEP Random           5 //                             HEP Random
  5 //                       --- Ranlux64Engine --      6 //                       --- Ranlux64Engine ---
  6 //                      class implementation f      7 //                      class implementation file
  7 // -------------------------------------------      8 // -----------------------------------------------------------------------
  8 // A double-precision implementation of the Ra      9 // A double-precision implementation of the RanluxEngine generator as 
  9 // decsribed by the notes of the original ranl     10 // decsribed by the notes of the original ranlux author (Martin Luscher)
 10 //                                                 11 //
 11 // See the note by Martin Luscher, December 19     12 // See the note by Martin Luscher, December 1997, entitiled
 12 // Double-precision implementation of the rand     13 // Double-precision implementation of the random number generator ranlux
 13 //                                                 14 //
 14 // ===========================================     15 // =======================================================================
 15 // Ken Smith      - Initial draft: 14th Jul 19     16 // Ken Smith      - Initial draft: 14th Jul 1998
 16 //                - Removed pow() from flat me     17 //                - Removed pow() from flat method 14th Jul 1998
 17 //                - Added conversion operators     18 //                - Added conversion operators:  6th Aug 1998
 18 //                                                 19 //
 19 // Mark Fischler  The following were modified      20 // Mark Fischler  The following were modified mostly to make the routine
 20 //      exactly match the Luscher algorithm in     21 //      exactly match the Luscher algorithm in generating 48-bit
 21 //      randoms:                                   22 //      randoms:
 22 // 9/9/98   - Substantial changes in what used     23 // 9/9/98   - Substantial changes in what used to be flat() to match
 23 //        algorithm in Luscher's ranlxd.c          24 //        algorithm in Luscher's ranlxd.c
 24 //      - Added update() method for 12 numbers     25 //      - Added update() method for 12 numbers, making flat() trivial
 25 //      - Added advance() method to hold the u     26 //      - Added advance() method to hold the unrolled loop for update
 26 //      - Distinction between three forms of s     27 //      - Distinction between three forms of seeding such that it
 27 //        is impossible to get same sequence f     28 //        is impossible to get same sequence from different forms -
 28 //        done by discarding some fraction of      29 //        done by discarding some fraction of one macro cycle which
 29 //        is different for the three cases         30 //        is different for the three cases
 30 //      - Change the misnomer "seed_table" to      31 //      - Change the misnomer "seed_table" to the more accurate 
 31 //        "randoms"                                32 //        "randoms"
 32 //      - Removed the no longer needed count12     33 //      - Removed the no longer needed count12, i_lag, j_lag, etc.
 33 //      - Corrected seed procedure which had b     34 //      - Corrected seed procedure which had been filling bits past
 34 //        2^-48.  This actually was very bad,      35 //        2^-48.  This actually was very bad, invalidating the
 35 //        number theory behind the proof that      36 //        number theory behind the proof that ranlxd is good.
 36 //      - Addition of 2**(-49) to generated nu     37 //      - Addition of 2**(-49) to generated number to prevent zero 
 37 //        from being returned; this does not a     38 //        from being returned; this does not affect the sequence 
 38 //        itself.                                  39 //        itself.
 39 //      - Corrected ecu seeding, which had bee     40 //      - Corrected ecu seeding, which had been supplying only 
 40 //        numbers less than 1/2.  This is prob     41 //        numbers less than 1/2.  This is probably moot.
 41 // 9/15/98    - Modified use of the various ex     42 // 9/15/98    - Modified use of the various exponents of 2
 42 //                  to avoid per-instance spac     43 //                  to avoid per-instance space overhead.  Note that these
 43 //        are initialized in setSeed, which EV     44 //        are initialized in setSeed, which EVERY constructor
 44 //        must invoke.                             45 //        must invoke.
 45 // J. Marraffino  - Remove dependence on hepSt     46 // J. Marraffino  - Remove dependence on hepString class  13 May 1999
 46 // M. Fischler    - In restore, checkFile for      47 // M. Fischler    - In restore, checkFile for file not found    03 Dec 2004
 47 // M. Fischler    - put get Methods for distri     48 // M. Fischler    - put get Methods for distrib instance save/restore 12/8/04    
 48 // M. Fischler    - split get() into tag valid     49 // M. Fischler    - split get() into tag validation and 
 49 //                  getState() for anonymous r     50 //                  getState() for anonymous restores           12/27/04    
 50 // M. Fischler    - put/get for vectors of ulo     51 // M. Fischler    - put/get for vectors of ulongs   3/14/05
 51 // M. Fischler    - State-saving using only in     52 // M. Fischler    - State-saving using only ints, for portability 4/12/05
 52 //                                                 53 //
 53 // ===========================================     54 // =======================================================================
 54                                                    55 
 55 #include "CLHEP/Random/Random.h"                   56 #include "CLHEP/Random/Random.h"
 56 #include "CLHEP/Random/Ranlux64Engine.h"           57 #include "CLHEP/Random/Ranlux64Engine.h"
 57 #include "CLHEP/Random/engineIDulong.h"            58 #include "CLHEP/Random/engineIDulong.h"
 58 #include "CLHEP/Random/DoubConv.h"                 59 #include "CLHEP/Random/DoubConv.h"
 59 #include "CLHEP/Utility/atomic_int.h"          <<  60 #include <string.h> // for strcmp
 60                                                << 
 61 #include <atomic>                              << 
 62 #include <cstdlib>  // for std::abs(int)           61 #include <cstdlib>  // for std::abs(int)
 63 #include <iostream>                            << 
 64 #include <limits>       // for numeric_limits      62 #include <limits>       // for numeric_limits
 65 #include <string.h> // for strcmp              << 
 66 #include <vector>                              << 
 67                                                    63 
 68 namespace CLHEP {                                  64 namespace CLHEP {
 69                                                    65 
 70 namespace {                                    <<  66 static const int MarkerLen = 64; // Enough room to hold a begin or end marker. 
 71   // Number of instances with automatic seed s << 
 72   CLHEP_ATOMIC_INT_TYPE numberOfEngines(0);    << 
 73                                                    67 
 74   // Maximum index into the seed table         << 
 75   const int maxIndex = 215;                    << 
 76 }                                              << 
 77                                                    68 
 78 static const int MarkerLen = 64; // Enough roo <<  69 // Number of instances with automatic seed selection
                                                   >>  70 int Ranlux64Engine::numEngines = 0;
 79                                                    71 
                                                   >>  72 // Maximum index into the seed table
                                                   >>  73 int Ranlux64Engine::maxIndex = 215;
 80                                                    74 
 81 #ifndef WIN32                                      75 #ifndef WIN32
 82 namespace detail {                                 76 namespace detail {
 83                                                    77 
 84 template< std::size_t n,                           78 template< std::size_t n,
 85           bool = n < std::size_t(std::numeric_     79           bool = n < std::size_t(std::numeric_limits<unsigned long>::digits) >
 86   struct do_right_shift;                           80   struct do_right_shift;
 87 template< std::size_t n >                          81 template< std::size_t n >
 88   struct do_right_shift<n,true>                    82   struct do_right_shift<n,true>
 89 {                                                  83 {
 90   unsigned long operator()(unsigned long value     84   unsigned long operator()(unsigned long value) { return value >> n; }
 91 };                                                 85 };
 92 template< std::size_t n >                          86 template< std::size_t n >
 93   struct do_right_shift<n,false>                   87   struct do_right_shift<n,false>
 94 {                                                  88 {
 95   unsigned long operator()(unsigned long) { re     89   unsigned long operator()(unsigned long) { return 0ul; }
 96 };                                                 90 };
 97                                                    91 
 98 template< std::size_t nbits >                      92 template< std::size_t nbits >
 99   unsigned long rshift( unsigned long value )      93   unsigned long rshift( unsigned long value )
100 { return do_right_shift<nbits>()(value); }         94 { return do_right_shift<nbits>()(value); }
101                                                    95 
102 } // namespace detail                              96 } // namespace detail
103 #endif                                             97 #endif
104                                                    98 
105 std::string Ranlux64Engine::name() const {retu     99 std::string Ranlux64Engine::name() const {return "Ranlux64Engine";}
106                                                   100 
107 Ranlux64Engine::Ranlux64Engine()                  101 Ranlux64Engine::Ranlux64Engine()
108 : HepRandomEngine()                               102 : HepRandomEngine()
109 {                                                 103 {
110    luxury = 1;                                    104    luxury = 1;
111    int numEngines = numberOfEngines++;         << 
112    int cycle    = std::abs(int(numEngines/maxI    105    int cycle    = std::abs(int(numEngines/maxIndex));
113    int curIndex = std::abs(int(numEngines%maxI    106    int curIndex = std::abs(int(numEngines%maxIndex));
114                                                << 107    numEngines +=1;
115    long mask = ((cycle & 0x007fffff) << 8);       108    long mask = ((cycle & 0x007fffff) << 8);
116    long seedlist[2];                              109    long seedlist[2];
117    HepRandom::getTheTableSeeds( seedlist, curI    110    HepRandom::getTheTableSeeds( seedlist, curIndex );
118    seedlist[0] ^= mask;                           111    seedlist[0] ^= mask;
119    seedlist[1] = 0;                               112    seedlist[1] = 0;
120                                                   113 
121    setSeeds(seedlist, luxury);                    114    setSeeds(seedlist, luxury);
122    advance ( 8 );     // Discard some iteratio    115    advance ( 8 );     // Discard some iterations and ensure that
123         // this sequence won't match one where    116         // this sequence won't match one where seeds 
124         // were provided.                         117         // were provided.
125 }                                                 118 }
126                                                   119 
127 Ranlux64Engine::Ranlux64Engine(long seed, int     120 Ranlux64Engine::Ranlux64Engine(long seed, int lux)
128 : HepRandomEngine()                               121 : HepRandomEngine()
129 {                                                 122 {
130    luxury = lux;                                  123    luxury = lux;
131    long seedlist[2]={seed,0};                     124    long seedlist[2]={seed,0};
132    setSeeds(seedlist, lux);                       125    setSeeds(seedlist, lux);
133    advance ( 2*lux + 1 );   // Discard some it    126    advance ( 2*lux + 1 );   // Discard some iterations to use a different 
134         // point in the sequence.                 127         // point in the sequence.  
135 }                                                 128 }
136                                                   129 
137 Ranlux64Engine::Ranlux64Engine(int rowIndex, i    130 Ranlux64Engine::Ranlux64Engine(int rowIndex, int, int lux)
138 : HepRandomEngine()                               131 : HepRandomEngine()
139 {                                                 132 {
140    luxury = lux;                                  133    luxury = lux;
141    int cycle = std::abs(int(rowIndex/maxIndex)    134    int cycle = std::abs(int(rowIndex/maxIndex));
142    int   row = std::abs(int(rowIndex%maxIndex)    135    int   row = std::abs(int(rowIndex%maxIndex));
143    long mask = (( cycle & 0x000007ff ) << 20 )    136    long mask = (( cycle & 0x000007ff ) << 20 );
144    long seedlist[2];                              137    long seedlist[2]; 
145    HepRandom::getTheTableSeeds( seedlist, row     138    HepRandom::getTheTableSeeds( seedlist, row );
146    seedlist[0] ^= mask;                           139    seedlist[0] ^= mask;
147    seedlist[1]= 0;                                140    seedlist[1]= 0;
148    setSeeds(seedlist, lux);                       141    setSeeds(seedlist, lux);
149 }                                                 142 }
150                                                   143 
151 Ranlux64Engine::Ranlux64Engine( std::istream&     144 Ranlux64Engine::Ranlux64Engine( std::istream& is )
152 : HepRandomEngine()                               145 : HepRandomEngine()
153 {                                                 146 {
154   is >> *this;                                    147   is >> *this;
155 }                                                 148 }
156                                                   149 
157 Ranlux64Engine::~Ranlux64Engine() {}              150 Ranlux64Engine::~Ranlux64Engine() {}
158                                                   151 
159 double Ranlux64Engine::flat() {                   152 double Ranlux64Engine::flat() {
160   // Luscher improves the speed by computing s    153   // Luscher improves the speed by computing several numbers in a shot,
161   // in a manner similar to that of the Tauswo    154   // in a manner similar to that of the Tausworth in DualRand or the Hurd
162   // engines.  Thus, the real work is done in     155   // engines.  Thus, the real work is done in update().  Here we merely ensure
163   // that zero, which the algorithm can produc    156   // that zero, which the algorithm can produce, is never returned by flat().
164                                                   157 
165   if (index <= 0) update();                       158   if (index <= 0) update();
166   return randoms[--index] + twoToMinus_49();      159   return randoms[--index] + twoToMinus_49();
167 }                                                 160 }
168                                                   161 
169 void Ranlux64Engine::update() {                   162 void Ranlux64Engine::update() {
170   // Update the stash of twelve random numbers    163   // Update the stash of twelve random numbers.  
171   // When this routione is entered, index is a    164   // When this routione is entered, index is always 0.  The randoms 
172   // contains the last 12 numbers in the seque    165   // contains the last 12 numbers in the sequents:  s[0] is x[a+11], 
173   // s[1] is x[a+10] ... and s[11] is x[a] for    166   // s[1] is x[a+10] ... and s[11] is x[a] for some a.  Carry contains
174   // the last carry value (c[a+11]).              167   // the last carry value (c[a+11]).
175   //                                              168   //
176   // The recursion relation (3) in Luscher's n    169   // The recursion relation (3) in Luscher's note says 
177   //   delta[n] = x[n-s] = x[n-r] -c[n-1] or f    170   //   delta[n] = x[n-s] = x[n-r] -c[n-1] or for n=a+12,
178   //   delta[a+12] = x[a+7] - x[a] -c[a+11] wh    171   //   delta[a+12] = x[a+7] - x[a] -c[a+11] where we use r=12, s=5 per eqn. (7)
179   // This reduces to                              172   // This reduces to 
180   // s[11] = s[4] - s[11] - carry.                173   // s[11] = s[4] - s[11] - carry.
181   // The next number similarly will be given b    174   // The next number similarly will be given by s[10] = s[3] - s[10] - carry,
182   // and so forth until s[0] is filled.           175   // and so forth until s[0] is filled.
183   //                                              176   // 
184   // However, we need to skip 397, 202 or 109     177   // However, we need to skip 397, 202 or 109 numbers - these are not divisible 
185   // by 12 - to "fare well in the spectral tes    178   // by 12 - to "fare well in the spectral test".  
186                                                   179 
187   advance(pDozens);                               180   advance(pDozens);
188                                                   181 
189   // Since we wish at the end to have the 12 l    182   // Since we wish at the end to have the 12 last numbers in the order of 
190   // s[11] first, till s[0] last, we will have    183   // s[11] first, till s[0] last, we will have to do 1, 10, or 1 iterations 
191   // and then re-arrange to place to get the o    184   // and then re-arrange to place to get the oldest one in s[11].
192   // Generically, this will imply re-arranging    185   // Generically, this will imply re-arranging the s array at the end,
193   // but we can treat the special case of endI    186   // but we can treat the special case of endIters = 1 separately for superior
194   // efficiency in the cases of levels 0 and 2    187   // efficiency in the cases of levels 0 and 2.
195                                                   188 
196   double  y1;                                  << 189   register double  y1;
197                                                   190 
198   if ( endIters == 1 ) {    // Luxury levels 0    191   if ( endIters == 1 ) {    // Luxury levels 0 and 2 will go here
199     y1 = randoms[ 4] - randoms[11] - carry;       192     y1 = randoms[ 4] - randoms[11] - carry;
200     if ( y1 < 0.0 ) {                             193     if ( y1 < 0.0 ) {
201       y1 += 1.0;                                  194       y1 += 1.0;      
202       carry = twoToMinus_48();                    195       carry = twoToMinus_48();
203     } else {                                      196     } else {
204       carry = 0.0;                                197       carry = 0.0;
205     }                                             198     }
206     randoms[11] = randoms[10];                    199     randoms[11] = randoms[10];  
207     randoms[10] = randoms[ 9];                    200     randoms[10] = randoms[ 9];  
208     randoms[ 9] = randoms[ 8];                    201     randoms[ 9] = randoms[ 8];  
209     randoms[ 8] = randoms[ 7];                    202     randoms[ 8] = randoms[ 7];  
210     randoms[ 7] = randoms[ 6];                    203     randoms[ 7] = randoms[ 6];  
211     randoms[ 6] = randoms[ 5];                    204     randoms[ 6] = randoms[ 5];  
212     randoms[ 5] = randoms[ 4];                    205     randoms[ 5] = randoms[ 4];  
213     randoms[ 4] = randoms[ 3];                    206     randoms[ 4] = randoms[ 3];  
214     randoms[ 3] = randoms[ 2];                    207     randoms[ 3] = randoms[ 2];  
215     randoms[ 2] = randoms[ 1];                    208     randoms[ 2] = randoms[ 1];  
216     randoms[ 1] = randoms[ 0];                    209     randoms[ 1] = randoms[ 0];  
217     randoms[ 0] = y1;                             210     randoms[ 0] = y1;
218                                                   211 
219   } else {                                        212   } else {
220                                                   213 
221     int m, nr, ns;                                214     int m, nr, ns;
222     for ( m = 0, nr = 11, ns = 4; m < endIters    215     for ( m = 0, nr = 11, ns = 4; m < endIters; ++m, --nr ) {
223       y1 = randoms [ns] - randoms[nr] - carry;    216       y1 = randoms [ns] - randoms[nr] - carry;
224       if ( y1 < 0.0 ) {                           217       if ( y1 < 0.0 ) {
225         y1 += 1.0;                                218         y1 += 1.0;
226         carry = twoToMinus_48();                  219         carry = twoToMinus_48();
227       } else {                                    220       } else {
228         carry = 0.0;                              221         carry = 0.0;
229       }                                           222       }
230       randoms[nr] = y1;                           223       randoms[nr] = y1;
231       --ns;                                       224       --ns;
232       if ( ns < 0 ) {                             225       if ( ns < 0 ) {
233         ns = 11;                                  226         ns = 11;
234       }                                           227       }
235     } // loop on m                                228     } // loop on m
236                                                   229 
237     double temp[12];                              230     double temp[12];
238     for (m=0; m<12; m++) {                        231     for (m=0; m<12; m++) {
239       temp[m]=randoms[m];                         232       temp[m]=randoms[m];
240     }                                             233     }
241                                                   234 
242     ns = 11 - endIters;                           235     ns = 11 - endIters;
243     for (m=11; m>=0; --m) {                       236     for (m=11; m>=0; --m) {
244       randoms[m] = temp[ns];                      237       randoms[m] = temp[ns];
245       --ns;                                       238       --ns;
246       if ( ns < 0 ) {                             239       if ( ns < 0 ) {
247         ns = 11;                                  240         ns = 11;
248       }                                           241       }
249     }                                             242     } 
250                                                   243 
251   }                                               244   }
252                                                   245 
253   // Now when we return, there are 12 fresh us    246   // Now when we return, there are 12 fresh usable numbers in s[11] ... s[0]
254                                                   247 
255   index = 12;                                  << 248   index = 11;
256                                                   249 
257 } // update()                                     250 } // update()
258                                                   251 
259 void Ranlux64Engine::advance(int dozens) {        252 void Ranlux64Engine::advance(int dozens) {
260                                                   253 
261   double  y1, y2, y3;                          << 254   register double  y1, y2, y3;
262   double  cValue = twoToMinus_48();            << 255   register double  cValue = twoToMinus_48();
263   double  zero = 0.0;                          << 256   register double  zero = 0.0;
264   double  one  = 1.0;                          << 257   register double  one  = 1.0;
265                                                   258 
266     // Technical note:  We use Luscher's trick    259     // Technical note:  We use Luscher's trick to only do the
267     // carry subtraction when we really have t    260     // carry subtraction when we really have to.  Like him, we use 
268     // three registers instead of two so that     261     // three registers instead of two so that we avoid sequences
269     // like storing y1 then immediately replac    262     // like storing y1 then immediately replacing its value:
270     // some architectures lose time when this     263     // some architectures lose time when this is done.
271                                                   264 
272       // Luscher's ranlxd.c fills the stash go    265       // Luscher's ranlxd.c fills the stash going
273     // upward.  We fill it downward to save a     266     // upward.  We fill it downward to save a bit of time in the
274     // flat() routine at no cost later.  This     267     // flat() routine at no cost later.  This means that while
275     // Luscher's ir is jr+5, our n-r is (n-s)-    268     // Luscher's ir is jr+5, our n-r is (n-s)-5.  (Note that
276     // though ranlxd.c initializes ir and jr t    269     // though ranlxd.c initializes ir and jr to 11 and 7, ir as
277     // used is 5 more than jr because update i    270     // used is 5 more than jr because update is entered after 
278     // incrementing ir.)                          271     // incrementing ir.)  
279     //                                            272     //
280                                                   273 
281     // I have CAREFULLY checked that the algor    274     // I have CAREFULLY checked that the algorithms do match
282     // in all details.                            275     // in all details.
283                                                   276 
284   int k;                                          277   int k;
285   for ( k = dozens; k > 0; --k ) {                278   for ( k = dozens; k > 0; --k ) {
286                                                   279 
287     y1 = randoms[ 4] - randoms[11] - carry;       280     y1 = randoms[ 4] - randoms[11] - carry;
288                                                   281 
289     y2 = randoms[ 3] - randoms[10];               282     y2 = randoms[ 3] - randoms[10];
290     if ( y1 < zero ) {                            283     if ( y1 < zero ) {
291       y1 += one;                                  284       y1 += one;      
292       y2 -= cValue;                               285       y2 -= cValue;
293     }                                             286     }
294     randoms[11] = y1;                             287     randoms[11] = y1;
295                                                   288 
296     y3 = randoms[ 2] - randoms[ 9];               289     y3 = randoms[ 2] - randoms[ 9];
297     if ( y2 < zero ) {                            290     if ( y2 < zero ) {
298       y2 += one;                                  291       y2 += one;      
299       y3 -= cValue;                               292       y3 -= cValue;
300     }                                             293     }
301     randoms[10] = y2;                             294     randoms[10] = y2;
302                                                   295 
303     y1 = randoms[ 1] - randoms[ 8];               296     y1 = randoms[ 1] - randoms[ 8];
304     if ( y3 < zero ) {                            297     if ( y3 < zero ) {
305       y3 += one;                                  298       y3 += one;      
306       y1 -= cValue;                               299       y1 -= cValue;
307     }                                             300     }
308     randoms[ 9] = y3;                             301     randoms[ 9] = y3;
309                                                   302 
310     y2 = randoms[ 0] - randoms[ 7];               303     y2 = randoms[ 0] - randoms[ 7];
311     if ( y1 < zero ) {                            304     if ( y1 < zero ) {
312       y1 += one;                                  305       y1 += one;      
313       y2 -= cValue;                               306       y2 -= cValue;
314     }                                             307     }
315     randoms[ 8] = y1;                             308     randoms[ 8] = y1;
316                                                   309 
317     y3 = randoms[11] - randoms[ 6];               310     y3 = randoms[11] - randoms[ 6];
318     if ( y2 < zero ) {                            311     if ( y2 < zero ) {
319       y2 += one;                                  312       y2 += one;      
320       y3 -= cValue;                               313       y3 -= cValue;
321     }                                             314     }
322     randoms[ 7] = y2;                             315     randoms[ 7] = y2;
323                                                   316 
324     y1 = randoms[10] - randoms[ 5];               317     y1 = randoms[10] - randoms[ 5];
325     if ( y3 < zero ) {                            318     if ( y3 < zero ) {
326       y3 += one;                                  319       y3 += one;      
327       y1 -= cValue;                               320       y1 -= cValue;
328     }                                             321     }
329     randoms[ 6] = y3;                             322     randoms[ 6] = y3;
330                                                   323 
331     y2 = randoms[ 9] - randoms[ 4];               324     y2 = randoms[ 9] - randoms[ 4];
332     if ( y1 < zero ) {                            325     if ( y1 < zero ) {
333       y1 += one;                                  326       y1 += one;      
334       y2 -= cValue;                               327       y2 -= cValue;
335     }                                             328     }
336     randoms[ 5] = y1;                             329     randoms[ 5] = y1;
337                                                   330 
338     y3 = randoms[ 8] - randoms[ 3];               331     y3 = randoms[ 8] - randoms[ 3];
339     if ( y2 < zero ) {                            332     if ( y2 < zero ) {
340       y2 += one;                                  333       y2 += one;      
341       y3 -= cValue;                               334       y3 -= cValue;
342     }                                             335     }
343     randoms[ 4] = y2;                             336     randoms[ 4] = y2;
344                                                   337 
345     y1 = randoms[ 7] - randoms[ 2];               338     y1 = randoms[ 7] - randoms[ 2];
346     if ( y3 < zero ) {                            339     if ( y3 < zero ) {
347       y3 += one;                                  340       y3 += one;      
348       y1 -= cValue;                               341       y1 -= cValue;
349     }                                             342     }
350     randoms[ 3] = y3;                             343     randoms[ 3] = y3;
351                                                   344 
352     y2 = randoms[ 6] - randoms[ 1];               345     y2 = randoms[ 6] - randoms[ 1];
353     if ( y1 < zero ) {                            346     if ( y1 < zero ) {
354       y1 += one;                                  347       y1 += one;      
355       y2 -= cValue;                               348       y2 -= cValue;
356     }                                             349     }
357     randoms[ 2] = y1;                             350     randoms[ 2] = y1;
358                                                   351 
359     y3 = randoms[ 5] - randoms[ 0];               352     y3 = randoms[ 5] - randoms[ 0];
360     if ( y2 < zero ) {                            353     if ( y2 < zero ) {
361       y2 += one;                                  354       y2 += one;      
362       y3 -= cValue;                               355       y3 -= cValue;
363     }                                             356     }
364     randoms[ 1] = y2;                             357     randoms[ 1] = y2;
365                                                   358 
366     if ( y3 < zero ) {                            359     if ( y3 < zero ) {
367       y3 += one;                                  360       y3 += one;      
368       carry = cValue;                             361       carry = cValue;
369     }                                             362     }
370     randoms[ 0] = y3;                             363     randoms[ 0] = y3; 
371                                                   364 
372   } // End of major k loop doing 12 numbers at    365   } // End of major k loop doing 12 numbers at each cycle
373                                                   366 
374 } // advance(dozens)                              367 } // advance(dozens)
375                                                   368 
376 void Ranlux64Engine::flatArray(const int size,    369 void Ranlux64Engine::flatArray(const int size, double* vect) {
377   for( int i=0; i < size; ++i ) {                 370   for( int i=0; i < size; ++i ) {
378     vect[i] = flat();                             371     vect[i] = flat(); 
379   }                                               372   }
380 }                                                 373 }
381                                                   374 
382 void Ranlux64Engine::setSeed(long seed, int lu    375 void Ranlux64Engine::setSeed(long seed, int lux) {
383                                                   376 
384 // The initialization is carried out using a M    377 // The initialization is carried out using a Multiplicative
385 // Congruential generator using formula consta    378 // Congruential generator using formula constants of L'Ecuyer
386 // as described in "A review of pseudorandom n    379 // as described in "A review of pseudorandom number generators"
387 // (Fred James) published in Computer Physics     380 // (Fred James) published in Computer Physics Communications 60 (1990)
388 // pages 329-344                                  381 // pages 329-344
389                                                   382 
390   const int ecuyer_a(53668);                      383   const int ecuyer_a(53668);
391   const int ecuyer_b(40014);                      384   const int ecuyer_b(40014);
392   const int ecuyer_c(12211);                      385   const int ecuyer_c(12211);
393   const int ecuyer_d(2147483563);                 386   const int ecuyer_d(2147483563);
394                                                   387 
395   const int lux_levels[3] = {109, 202, 397};      388   const int lux_levels[3] = {109, 202, 397};
396   theSeed = seed;                                 389   theSeed = seed;
397                                                   390 
398   if( (lux > 2)||(lux < 0) ){                     391   if( (lux > 2)||(lux < 0) ){
399      pDiscard = (lux >= 12) ? (lux-12) : lux_l    392      pDiscard = (lux >= 12) ? (lux-12) : lux_levels[1];
400   }else{                                          393   }else{
401      pDiscard = lux_levels[luxury];               394      pDiscard = lux_levels[luxury];
402   }                                               395   }
403   pDozens  = pDiscard / 12;                       396   pDozens  = pDiscard / 12;
404   endIters = pDiscard % 12;                       397   endIters = pDiscard % 12;
405                                                   398 
406   long init_table[24];                            399   long init_table[24];
407   long next_seed = seed;                          400   long next_seed = seed;
408   long k_multiple;                                401   long k_multiple;
409   int i;                                          402   int i;
410   next_seed &= 0xffffffff;                        403   next_seed &= 0xffffffff;
411   while( next_seed >= ecuyer_d ) {                404   while( next_seed >= ecuyer_d ) {
412      next_seed -= ecuyer_d;                       405      next_seed -= ecuyer_d;
413   }                                               406   }
414                                                   407   
415   for(i = 0;i != 24;i++){                         408   for(i = 0;i != 24;i++){
416      k_multiple = next_seed / ecuyer_a;           409      k_multiple = next_seed / ecuyer_a;
417      next_seed = ecuyer_b * (next_seed - k_mul    410      next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a)
418                                        - k_mul    411                                        - k_multiple * ecuyer_c;
419      if(next_seed < 0) {                          412      if(next_seed < 0) {
420   next_seed += ecuyer_d;                          413   next_seed += ecuyer_d;
421      }                                            414      }
422      next_seed &= 0xffffffff;                     415      next_seed &= 0xffffffff;
423      init_table[i] = next_seed;                   416      init_table[i] = next_seed;
424   }                                               417   } 
425   // are we on a 64bit machine?                   418   // are we on a 64bit machine?
426   if( sizeof(long) >= 8 ) {                       419   if( sizeof(long) >= 8 ) {
427      int64_t topbits1, topbits2;               << 420      long topbits1, topbits2;
428 #ifdef WIN32                                      421 #ifdef WIN32
429      topbits1 = ( (int64_t) seed >> 32) & 0xff << 422      topbits1 = ( seed >> 32) & 0xffff ;
430      topbits2 = ( (int64_t) seed >> 48) & 0xff << 423      topbits2 = ( seed >> 48) & 0xffff ;
431 #else                                             424 #else
432      topbits1 = detail::rshift<32>(seed) & 0xf    425      topbits1 = detail::rshift<32>(seed) & 0xffff ;
433      topbits2 = detail::rshift<48>(seed) & 0xf    426      topbits2 = detail::rshift<48>(seed) & 0xffff ;
434 #endif                                            427 #endif
435      init_table[0] ^= topbits1;                   428      init_table[0] ^= topbits1;
436      init_table[2] ^= topbits2;                   429      init_table[2] ^= topbits2;
437      //std::cout << " init_table[0] " << init_    430      //std::cout << " init_table[0] " << init_table[0] << " from " << topbits1 << std::endl;
438      //std::cout << " init_table[2] " << init_    431      //std::cout << " init_table[2] " << init_table[2] << " from " << topbits2 << std::endl;
439   }                                               432   }   
440                                                   433 
441   for(i = 0;i < 12; i++){                         434   for(i = 0;i < 12; i++){
442      randoms[i] = (init_table[2*i  ]      ) *     435      randoms[i] = (init_table[2*i  ]      ) * 2.0 * twoToMinus_32() +
443                   (init_table[2*i+1] >> 15) *     436                   (init_table[2*i+1] >> 15) * twoToMinus_48();
444      //if( randoms[i] < 0. || randoms[i]  > 1.    437      //if( randoms[i] < 0. || randoms[i]  > 1. ) {
445      //std::cout << "setSeed:  init_table " <<    438      //std::cout << "setSeed:  init_table " << init_table[2*i  ] << std::endl;
446      //std::cout << "setSeed:  init_table " <<    439      //std::cout << "setSeed:  init_table " << init_table[2*i+1] << std::endl;
447      //std::cout << "setSeed:  random " << i <    440      //std::cout << "setSeed:  random " << i << " is " << randoms[i] << std::endl;
448      //}                                          441      //}
449   }                                               442   }
450                                                   443 
451   carry = 0.0;                                    444   carry = 0.0;
452   if ( randoms[11] == 0. ) carry = twoToMinus_    445   if ( randoms[11] == 0. ) carry = twoToMinus_48();
453   // Perform an update before returning the fi << 446   index = 11;
454   index = -1;                                  << 
455                                                   447 
456 } // setSeed()                                    448 } // setSeed()
457                                                   449 
458 void Ranlux64Engine::setSeeds(const long * see    450 void Ranlux64Engine::setSeeds(const long * seeds, int lux) {
459 // old code only uses the first long in seeds     451 // old code only uses the first long in seeds
460 //  setSeed( *seeds ? *seeds : 32767, lux );      452 //  setSeed( *seeds ? *seeds : 32767, lux );
461 //  theSeeds = seeds;                             453 //  theSeeds = seeds;
462                                                   454 
463 // using code from Ranlux - even those are 32b    455 // using code from Ranlux - even those are 32bit seeds, 
464 // that is good enough to completely different    456 // that is good enough to completely differentiate the sequences
465                                                   457 
466    const int ecuyer_a = 53668;                    458    const int ecuyer_a = 53668;
467    const int ecuyer_b = 40014;                    459    const int ecuyer_b = 40014;
468    const int ecuyer_c = 12211;                    460    const int ecuyer_c = 12211;
469    const int ecuyer_d = 2147483563;               461    const int ecuyer_d = 2147483563;
470                                                   462 
471    const int lux_levels[3] = {109, 202, 397};     463    const int lux_levels[3] = {109, 202, 397};
472    const long *seedptr;                           464    const long *seedptr; 
473                                                   465 
474    theSeeds = seeds;                              466    theSeeds = seeds;
475    seedptr  = seeds;                              467    seedptr  = seeds;
476                                                   468  
477    if(seeds == 0){                                469    if(seeds == 0){
478       setSeed(theSeed,lux);                       470       setSeed(theSeed,lux);
479       theSeeds = &theSeed;                        471       theSeeds = &theSeed;
480       return;                                     472       return;
481    }                                              473    }
482                                                   474 
483    theSeed = *seeds;                              475    theSeed = *seeds;
484                                                   476 
485 // number of additional random numbers that ne    477 // number of additional random numbers that need to be 'thrown away'
486 // every 24 numbers is set using luxury level     478 // every 24 numbers is set using luxury level variable.
487                                                   479 
488   if( (lux > 2)||(lux < 0) ){                     480   if( (lux > 2)||(lux < 0) ){
489      pDiscard = (lux >= 12) ? (lux-12) : lux_l    481      pDiscard = (lux >= 12) ? (lux-12) : lux_levels[1];
490   }else{                                          482   }else{
491      pDiscard = lux_levels[luxury];               483      pDiscard = lux_levels[luxury];
492   }                                               484   }
493   pDozens  = pDiscard / 12;                       485   pDozens  = pDiscard / 12;
494   endIters = pDiscard % 12;                       486   endIters = pDiscard % 12;
495                                                   487 
496   long init_table[24];                            488   long init_table[24];
497   long next_seed = *seeds;                        489   long next_seed = *seeds;
498   long k_multiple;                                490   long k_multiple;
499   int i;                                          491   int i;
500                                                   492       
501   for( i = 0;(i != 24)&&(*seedptr != 0);i++){     493   for( i = 0;(i != 24)&&(*seedptr != 0);i++){
502       init_table[i] =  *seedptr & 0xffffffff;     494       init_table[i] =  *seedptr & 0xffffffff;
503       seedptr++;                                  495       seedptr++;
504   }                                               496   }          
505                                                   497 
506   if(i != 24){                                    498   if(i != 24){
507      next_seed = init_table[i-1];                 499      next_seed = init_table[i-1];
508      for(;i != 24;i++){                           500      for(;i != 24;i++){
509   k_multiple = next_seed / ecuyer_a;              501   k_multiple = next_seed / ecuyer_a;
510   next_seed = ecuyer_b * (next_seed - k_multip    502   next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a)
511                                     - k_multip    503                                     - k_multiple * ecuyer_c;
512   if(next_seed < 0) {                             504   if(next_seed < 0) {
513      next_seed += ecuyer_d;                       505      next_seed += ecuyer_d;
514   }                                               506   }
515   next_seed &= 0xffffffff;                        507   next_seed &= 0xffffffff;
516   init_table[i] = next_seed;                      508   init_table[i] = next_seed;
517      }                                            509      }    
518   }                                               510   }
519                                                   511 
520   for(i = 0;i < 12; i++){                         512   for(i = 0;i < 12; i++){
521      randoms[i] = (init_table[2*i  ]      ) *     513      randoms[i] = (init_table[2*i  ]      ) * 2.0 * twoToMinus_32() +
522                   (init_table[2*i+1] >> 15) *     514                   (init_table[2*i+1] >> 15) * twoToMinus_48();
523   }                                               515   }
524                                                   516 
525   carry = 0.0;                                    517   carry = 0.0;
526   if ( randoms[11] == 0. ) carry = twoToMinus_    518   if ( randoms[11] == 0. ) carry = twoToMinus_48();
527   // Perform an update before returning the fi << 519   index = 11;
528   index = -1;                                  << 
529                                                   520 
530 }                                                 521 }
531                                                   522 
532 void Ranlux64Engine::saveStatus( const char fi    523 void Ranlux64Engine::saveStatus( const char filename[] ) const
533 {                                                 524 {
534    std::ofstream outFile( filename, std::ios::    525    std::ofstream outFile( filename, std::ios::out ) ;
535   if (!outFile.bad()) {                           526   if (!outFile.bad()) {
536     outFile << "Uvec\n";                          527     outFile << "Uvec\n";
537     std::vector<unsigned long> v = put();         528     std::vector<unsigned long> v = put();
538     for (unsigned int i=0; i<v.size(); ++i) {     529     for (unsigned int i=0; i<v.size(); ++i) {
539       outFile << v[i] << "\n";                    530       outFile << v[i] << "\n";
540     }                                             531     }
541   }                                               532   }
542 }                                                 533 }
543                                                   534 
544 void Ranlux64Engine::restoreStatus( const char    535 void Ranlux64Engine::restoreStatus( const char filename[] )
545 {                                                 536 {
546    std::ifstream inFile( filename, std::ios::i    537    std::ifstream inFile( filename, std::ios::in);
547    if (!checkFile ( inFile, filename, engineNa    538    if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) {
548      std::cerr << "  -- Engine state remains u    539      std::cerr << "  -- Engine state remains unchanged\n";
549      return;                                      540      return;
550    }                                              541    }
551   if ( possibleKeywordInput ( inFile, "Uvec",     542   if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) {
552     std::vector<unsigned long> v;                 543     std::vector<unsigned long> v;
553     unsigned long xin;                            544     unsigned long xin;
554     for (unsigned int ivec=0; ivec < VECTOR_ST    545     for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
555       inFile >> xin;                              546       inFile >> xin;
556       if (!inFile) {                              547       if (!inFile) {
557         inFile.clear(std::ios::badbit | inFile    548         inFile.clear(std::ios::badbit | inFile.rdstate());
558         std::cerr << "\nJamesRandom state (vec    549         std::cerr << "\nJamesRandom state (vector) description improper."
559          << "\nrestoreStatus has failed."         550          << "\nrestoreStatus has failed."
560          << "\nInput stream is probably mispos    551          << "\nInput stream is probably mispositioned now." << std::endl;
561         return;                                   552         return;
562       }                                           553       }
563       v.push_back(xin);                           554       v.push_back(xin);
564     }                                             555     }
565     getState(v);                                  556     getState(v);
566     return;                                       557     return;
567   }                                               558   }
568                                                   559 
569    if (!inFile.bad() && !inFile.eof()) {          560    if (!inFile.bad() && !inFile.eof()) {
570 //     inFile >> theSeed;  removed -- encompas    561 //     inFile >> theSeed;  removed -- encompased by possibleKeywordInput
571      for (int i=0; i<12; ++i) {                   562      for (int i=0; i<12; ++i) {
572        inFile >> randoms[i];                      563        inFile >> randoms[i];
573      }                                            564      }
574      inFile >> carry; inFile >> index;            565      inFile >> carry; inFile >> index;
575      inFile >> luxury; inFile >> pDiscard;        566      inFile >> luxury; inFile >> pDiscard;
576      pDozens  = pDiscard / 12;                    567      pDozens  = pDiscard / 12;
577      endIters = pDiscard % 12;                    568      endIters = pDiscard % 12;
578    }                                              569    }
579 }                                                 570 }
580                                                   571 
581 void Ranlux64Engine::showStatus() const           572 void Ranlux64Engine::showStatus() const
582 {                                                 573 {
583    std::cout << std::endl;                        574    std::cout << std::endl;
584    std::cout << "--------- Ranlux engine statu    575    std::cout << "--------- Ranlux engine status ---------" << std::endl;
585    std::cout << " Initial seed = " << theSeed     576    std::cout << " Initial seed = " << theSeed << std::endl;
586    std::cout << " randoms[] = ";                  577    std::cout << " randoms[] = ";
587    for (int i=0; i<12; ++i) {                     578    for (int i=0; i<12; ++i) {
588      std::cout << randoms[i] << std::endl;        579      std::cout << randoms[i] << std::endl;
589    }                                              580    }
590    std::cout << std::endl;                        581    std::cout << std::endl;
591    std::cout << " carry = " << carry << ", ind    582    std::cout << " carry = " << carry << ", index = " << index << std::endl;
592    std::cout << " luxury = " << luxury << " pD    583    std::cout << " luxury = " << luxury << " pDiscard = " 
593             << pDiscard << std::endl;             584             << pDiscard << std::endl;
594    std::cout << "-----------------------------    585    std::cout << "----------------------------------------" << std::endl;
595 }                                                 586 }
596                                                   587 
597 std::ostream & Ranlux64Engine::put( std::ostre    588 std::ostream & Ranlux64Engine::put( std::ostream& os ) const
598 {                                                 589 {
599    char beginMarker[] = "Ranlux64Engine-begin"    590    char beginMarker[] = "Ranlux64Engine-begin";
600   os << beginMarker << "\nUvec\n";                591   os << beginMarker << "\nUvec\n";
601   std::vector<unsigned long> v = put();           592   std::vector<unsigned long> v = put();
602   for (unsigned int i=0; i<v.size(); ++i) {       593   for (unsigned int i=0; i<v.size(); ++i) {
603      os <<  v[i] <<  "\n";                        594      os <<  v[i] <<  "\n";
604   }                                               595   }
605   return os;                                      596   return os;  
606 }                                                 597 }
607                                                   598 
608 std::vector<unsigned long> Ranlux64Engine::put    599 std::vector<unsigned long> Ranlux64Engine::put () const {
609   std::vector<unsigned long> v;                   600   std::vector<unsigned long> v;
610   v.push_back (engineIDulong<Ranlux64Engine>()    601   v.push_back (engineIDulong<Ranlux64Engine>());
611   std::vector<unsigned long> t;                   602   std::vector<unsigned long> t;
612   for (int i=0; i<12; ++i) {                      603   for (int i=0; i<12; ++i) {
613     t = DoubConv::dto2longs(randoms[i]);          604     t = DoubConv::dto2longs(randoms[i]);
614     v.push_back(t[0]); v.push_back(t[1]);         605     v.push_back(t[0]); v.push_back(t[1]);
615   }                                               606   }
616   t = DoubConv::dto2longs(carry);                 607   t = DoubConv::dto2longs(carry);
617   v.push_back(t[0]); v.push_back(t[1]);           608   v.push_back(t[0]); v.push_back(t[1]);
618   v.push_back(static_cast<unsigned long>(index    609   v.push_back(static_cast<unsigned long>(index));
619   v.push_back(static_cast<unsigned long>(luxur    610   v.push_back(static_cast<unsigned long>(luxury));
620   v.push_back(static_cast<unsigned long>(pDisc    611   v.push_back(static_cast<unsigned long>(pDiscard));
621   return v;                                       612   return v;
622 }                                                 613 }
623                                                   614 
624 std::istream & Ranlux64Engine::get ( std::istr    615 std::istream & Ranlux64Engine::get ( std::istream& is )
625 {                                                 616 {
626   char beginMarker [MarkerLen];                   617   char beginMarker [MarkerLen];
627   is >> std::ws;                                  618   is >> std::ws;
628   is.width(MarkerLen);  // causes the next rea    619   is.width(MarkerLen);  // causes the next read to the char* to be <=
629       // that many bytes, INCLUDING A TERMINAT    620       // that many bytes, INCLUDING A TERMINATION \0 
630       // (Stroustrup, section 21.3.2)             621       // (Stroustrup, section 21.3.2)
631   is >> beginMarker;                              622   is >> beginMarker;
632   if (strcmp(beginMarker,"Ranlux64Engine-begin    623   if (strcmp(beginMarker,"Ranlux64Engine-begin")) {
633      is.clear(std::ios::badbit | is.rdstate())    624      is.clear(std::ios::badbit | is.rdstate());
634      std::cerr << "\nInput stream mispositione    625      std::cerr << "\nInput stream mispositioned or"
635          << "\nRanlux64Engine state descriptio    626          << "\nRanlux64Engine state description missing or"
636          << "\nwrong engine type found." << st    627          << "\nwrong engine type found." << std::endl;
637      return is;                                   628      return is;
638   }                                               629   }
639   return getState(is);                            630   return getState(is);
640 }                                                 631 }
641                                                   632 
642 std::string Ranlux64Engine::beginTag ( )  {       633 std::string Ranlux64Engine::beginTag ( )  { 
643   return "Ranlux64Engine-begin";                  634   return "Ranlux64Engine-begin"; 
644 }                                                 635 }
645                                                   636 
646 std::istream & Ranlux64Engine::getState ( std:    637 std::istream & Ranlux64Engine::getState ( std::istream& is )
647 {                                                 638 {
648   if ( possibleKeywordInput ( is, "Uvec", theS    639   if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) {
649     std::vector<unsigned long> v;                 640     std::vector<unsigned long> v;
650     unsigned long uu;                             641     unsigned long uu;
651     for (unsigned int ivec=0; ivec < VECTOR_ST    642     for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
652       is >> uu;                                   643       is >> uu;
653       if (!is) {                                  644       if (!is) {
654         is.clear(std::ios::badbit | is.rdstate    645         is.clear(std::ios::badbit | is.rdstate());
655         std::cerr << "\nRanlux64Engine state (    646         std::cerr << "\nRanlux64Engine state (vector) description improper."
656     << "\ngetState() has failed."                 647     << "\ngetState() has failed."
657          << "\nInput stream is probably mispos    648          << "\nInput stream is probably mispositioned now." << std::endl;
658         return is;                                649         return is;
659       }                                           650       }
660       v.push_back(uu);                            651       v.push_back(uu);
661     }                                             652     }
662     getState(v);                                  653     getState(v);
663     return (is);                                  654     return (is);
664   }                                               655   }
665                                                   656 
666 //  is >> theSeed;  Removed, encompassed by po    657 //  is >> theSeed;  Removed, encompassed by possibleKeywordInput()
667                                                   658 
668   char endMarker   [MarkerLen];                   659   char endMarker   [MarkerLen];
669   for (int i=0; i<12; ++i) {                      660   for (int i=0; i<12; ++i) {
670      is >> randoms[i];                            661      is >> randoms[i];
671   }                                               662   }
672   is >> carry; is >> index;                       663   is >> carry; is >> index;
673   is >> luxury; is >> pDiscard;                   664   is >> luxury; is >> pDiscard;
674   pDozens  = pDiscard / 12;                       665   pDozens  = pDiscard / 12;
675   endIters = pDiscard % 12;                       666   endIters = pDiscard % 12;
676   is >> std::ws;                                  667   is >> std::ws;
677   is.width(MarkerLen);                            668   is.width(MarkerLen);  
678   is >> endMarker;                                669   is >> endMarker;
679   if (strcmp(endMarker,"Ranlux64Engine-end"))     670   if (strcmp(endMarker,"Ranlux64Engine-end")) {
680      is.clear(std::ios::badbit | is.rdstate())    671      is.clear(std::ios::badbit | is.rdstate());
681      std::cerr << "\nRanlux64Engine state desc    672      std::cerr << "\nRanlux64Engine state description incomplete."
682          << "\nInput stream is probably mispos    673          << "\nInput stream is probably mispositioned now." << std::endl;
683      return is;                                   674      return is;
684   }                                               675   }
685   return is;                                      676   return is;
686 }                                                 677 }
687                                                   678 
688 bool Ranlux64Engine::get (const std::vector<un    679 bool Ranlux64Engine::get (const std::vector<unsigned long> & v) {
689   if ((v[0] & 0xffffffffUL) != engineIDulong<R    680   if ((v[0] & 0xffffffffUL) != engineIDulong<Ranlux64Engine>()) {
690     std::cerr <<                                  681     std::cerr << 
691       "\nRanlux64Engine get:state vector has w    682       "\nRanlux64Engine get:state vector has wrong ID word - state unchanged\n";
692     return false;                                 683     return false;
693   }                                               684   }
694   return getState(v);                             685   return getState(v);
695 }                                                 686 }
696                                                   687 
697 bool Ranlux64Engine::getState (const std::vect    688 bool Ranlux64Engine::getState (const std::vector<unsigned long> & v) {
698   if (v.size() != VECTOR_STATE_SIZE ) {           689   if (v.size() != VECTOR_STATE_SIZE ) {
699     std::cerr <<                                  690     std::cerr << 
700       "\nRanlux64Engine get:state vector has w    691       "\nRanlux64Engine get:state vector has wrong length - state unchanged\n";
701     return false;                                 692     return false;
702   }                                               693   }
703   std::vector<unsigned long> t(2);                694   std::vector<unsigned long> t(2);
704   for (int i=0; i<12; ++i) {                      695   for (int i=0; i<12; ++i) {
705     t[0] = v[2*i+1]; t[1] = v[2*i+2];             696     t[0] = v[2*i+1]; t[1] = v[2*i+2];
706     randoms[i] = DoubConv::longs2double(t);       697     randoms[i] = DoubConv::longs2double(t);
707   }                                               698   }
708   t[0] = v[25]; t[1] = v[26];                     699   t[0] = v[25]; t[1] = v[26];
709   carry    = DoubConv::longs2double(t);           700   carry    = DoubConv::longs2double(t);
710   index    = (int)v[27];                       << 701   index    = v[27];
711   luxury   = (int)v[28];                       << 702   luxury   = v[28];
712   pDiscard = (int)v[29];                       << 703   pDiscard = v[29]; 
713   return true;                                    704   return true;
714 }                                                 705 }
715                                                   706 
716 }  // namespace CLHEP                             707 }  // namespace CLHEP
717                                                   708