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

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