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