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