1 // $Id: RanluxEngine.cc,v 1.6 2010/06/16 17:24:53 garren Exp $
2 // -*- C++ -*-
3 //
4 // -----------------------------------------------------------------------
5 // HEP Random
6 // --- RanluxEngine ---
7 // class implementation file
8 // -----------------------------------------------------------------------
9 // This file is part of Geant4 (simulation toolkit for HEP).
10 //
11 // Ranlux random number generator originally implemented in FORTRAN77
12 // by Fred James as part of the MATHLIB HEP library.
13 // 'RanluxEngine' is designed to fit into the CLHEP random number
14 // class structure.
15
16 // ===============================================================
17 // Adeyemi Adesanya - Created: 6th November 1995
18 // Gabriele Cosmo - Adapted & Revised: 22nd November 1995
19 // Adeyemi Adesanya - Added setSeeds() method: 2nd February 1996
20 // Gabriele Cosmo - Added flatArray() method: 8th February 1996
21 // - Minor corrections: 31st October 1996
22 // - Added methods for engine status: 19th November 1996
23 // - Fixed bug in setSeeds(): 15th September 1997
24 // J.Marraffino - Added stream operators and related constructor.
25 // Added automatic seed selection from seed table and
26 // engine counter: 14th Feb 1998
27 // - Fixed bug: setSeeds() requires a zero terminated
28 // array of seeds: 19th Feb 1998
29 // Ken Smith - Added conversion operators: 6th Aug 1998
30 // J. Marraffino - Remove dependence on hepString class 13 May 1999
31 // M. Fischler - In restore, checkFile for file not found 03 Dec 2004
32 // M. Fischler - Methods put, getfor instance save/restore 12/8/04
33 // M. Fischler - split get() into tag validation and
34 // getState() for anonymous restores 12/27/04
35 // M. Fischler - put/get for vectors of ulongs 3/14/05
36 // M. Fischler - State-saving using only ints, for portability 4/12/05
37 //
38 // ===============================================================
39
40 #include "CLHEP/Random/defs.h"
41 #include "CLHEP/Random/Random.h"
42 #include "CLHEP/Random/RanluxEngine.h"
43 #include "CLHEP/Random/engineIDulong.h"
44 #include "CLHEP/Utility/atomic_int.h"
45
46 #include <string.h> // for strcmp
47 #include <cstdlib> // for std::abs(int)
48
49 #ifdef TRACE_IO
50 #include "CLHEP/Random/DoubConv.hh"
51 bool flat_trace = false;
52 #endif
53
54 namespace CLHEP {
55
56 namespace {
57 // Number of instances with automatic seed selection
58 CLHEP_ATOMIC_INT_TYPE numberOfEngines(0);
59
60 // Maximum index into the seed table
61 const int maxIndex = 215;
62 }
63
64 static const int MarkerLen = 64; // Enough room to hold a begin or end marker.
65
name() const66 std::string RanluxEngine::name() const {return "RanluxEngine";}
67
RanluxEngine(long seed,int lux)68 RanluxEngine::RanluxEngine(long seed, int lux)
69 : HepRandomEngine()
70 {
71 long seedlist[2]={0,0};
72
73 luxury = lux;
74 setSeed(seed, luxury);
75
76 // setSeeds() wants a zero terminated array!
77 seedlist[0]=theSeed;
78 seedlist[1]=0;
79 setSeeds(seedlist, luxury);
80 }
81
RanluxEngine()82 RanluxEngine::RanluxEngine()
83 : HepRandomEngine()
84 {
85 long seed;
86 long seedlist[2]={0,0};
87
88 luxury = 3;
89 int numEngines = numberOfEngines++;
90 int cycle = std::abs(int(numEngines/maxIndex));
91 int curIndex = std::abs(int(numEngines%maxIndex));
92
93 long mask = ((cycle & 0x007fffff) << 8);
94 HepRandom::getTheTableSeeds( seedlist, curIndex );
95 seed = seedlist[0]^mask;
96 setSeed(seed, luxury);
97
98 // setSeeds() wants a zero terminated array!
99 seedlist[0]=theSeed;
100 seedlist[1]=0;
101 setSeeds(seedlist, luxury);
102 }
103
RanluxEngine(int rowIndex,int colIndex,int lux)104 RanluxEngine::RanluxEngine(int rowIndex, int colIndex, int lux)
105 : HepRandomEngine()
106 {
107 long seed;
108 long seedlist[2]={0,0};
109
110 luxury = lux;
111 int cycle = std::abs(int(rowIndex/maxIndex));
112 int row = std::abs(int(rowIndex%maxIndex));
113 int col = std::abs(int(colIndex%2));
114 long mask = (( cycle & 0x000007ff ) << 20 );
115 HepRandom::getTheTableSeeds( seedlist, row );
116 seed = ( seedlist[col] )^mask;
117 setSeed(seed, luxury);
118
119 // setSeeds() wants a zero terminated array!
120 seedlist[0]=theSeed;
121 seedlist[1]=0;
122 setSeeds(seedlist, luxury);
123 }
124
RanluxEngine(std::istream & is)125 RanluxEngine::RanluxEngine( std::istream& is )
126 : HepRandomEngine()
127 {
128 is >> *this;
129 }
130
~RanluxEngine()131 RanluxEngine::~RanluxEngine() {}
132
setSeed(long seed,int lux)133 void RanluxEngine::setSeed(long seed, int lux) {
134
135 // The initialisation is carried out using a Multiplicative
136 // Congruential generator using formula constants of L'Ecuyer
137 // as described in "A review of pseudorandom number generators"
138 // (Fred James) published in Computer Physics Communications 60 (1990)
139 // pages 329-344
140
141 const int ecuyer_a = 53668;
142 const int ecuyer_b = 40014;
143 const int ecuyer_c = 12211;
144 const int ecuyer_d = 2147483563;
145
146 const int lux_levels[5] = {0,24,73,199,365};
147
148 long int_seed_table[24];
149 long next_seed = seed;
150 long k_multiple;
151 int i;
152
153 // number of additional random numbers that need to be 'thrown away'
154 // every 24 numbers is set using luxury level variable.
155
156 theSeed = seed;
157 if( (lux > 4)||(lux < 0) ){
158 if(lux >= 24){
159 nskip = lux - 24;
160 }else{
161 nskip = lux_levels[3]; // corresponds to default luxury level
162 }
163 }else{
164 luxury = lux;
165 nskip = lux_levels[luxury];
166 }
167
168
169 for(i = 0;i != 24;i++){
170 k_multiple = next_seed / ecuyer_a;
171 next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a)
172 - k_multiple * ecuyer_c ;
173 if(next_seed < 0)next_seed += ecuyer_d;
174 int_seed_table[i] = next_seed % int_modulus;
175 }
176
177 for(i = 0;i != 24;i++)
178 float_seed_table[i] = int_seed_table[i] * mantissa_bit_24();
179
180 i_lag = 23;
181 j_lag = 9;
182 carry = 0. ;
183
184 if( float_seed_table[23] == 0. ) carry = mantissa_bit_24();
185
186 count24 = 0;
187 }
188
setSeeds(const long * seeds,int lux)189 void RanluxEngine::setSeeds(const long *seeds, int lux) {
190
191 const int ecuyer_a = 53668;
192 const int ecuyer_b = 40014;
193 const int ecuyer_c = 12211;
194 const int ecuyer_d = 2147483563;
195
196 const int lux_levels[5] = {0,24,73,199,365};
197 int i;
198 long int_seed_table[24];
199 long k_multiple,next_seed;
200 const long *seedptr;
201
202 theSeeds = seeds;
203 seedptr = seeds;
204
205 if(seeds == 0){
206 setSeed(theSeed,lux);
207 theSeeds = &theSeed;
208 return;
209 }
210
211 theSeed = *seeds;
212
213 // number of additional random numbers that need to be 'thrown away'
214 // every 24 numbers is set using luxury level variable.
215
216 if( (lux > 4)||(lux < 0) ){
217 if(lux >= 24){
218 nskip = lux - 24;
219 }else{
220 nskip = lux_levels[3]; // corresponds to default luxury level
221 }
222 }else{
223 luxury = lux;
224 nskip = lux_levels[luxury];
225 }
226
227 for( i = 0;(i != 24)&&(*seedptr != 0);i++){
228 int_seed_table[i] = *seedptr % int_modulus;
229 seedptr++;
230 }
231
232 if(i != 24){
233 next_seed = int_seed_table[i-1];
234 for(;i != 24;i++){
235 k_multiple = next_seed / ecuyer_a;
236 next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a)
237 - k_multiple * ecuyer_c ;
238 if(next_seed < 0)next_seed += ecuyer_d;
239 int_seed_table[i] = next_seed % int_modulus;
240 }
241 }
242
243 for(i = 0;i != 24;i++)
244 float_seed_table[i] = int_seed_table[i] * mantissa_bit_24();
245
246 i_lag = 23;
247 j_lag = 9;
248 carry = 0. ;
249
250 if( float_seed_table[23] == 0. ) carry = mantissa_bit_24();
251
252 count24 = 0;
253 }
254
saveStatus(const char filename[]) const255 void RanluxEngine::saveStatus( const char filename[] ) const
256 {
257 std::ofstream outFile( filename, std::ios::out ) ;
258 if (!outFile.bad()) {
259 outFile << "Uvec\n";
260 std::vector<unsigned long> v = put();
261 #ifdef TRACE_IO
262 std::cout << "Result of v = put() is:\n";
263 #endif
264 for (unsigned int i=0; i<v.size(); ++i) {
265 outFile << v[i] << "\n";
266 #ifdef TRACE_IO
267 std::cout << v[i] << " ";
268 if (i%6==0) std::cout << "\n";
269 #endif
270 }
271 #ifdef TRACE_IO
272 std::cout << "\n";
273 #endif
274 }
275 #ifdef REMOVED
276 if (!outFile.bad()) {
277 outFile << theSeed << std::endl;
278 for (int i=0; i<24; ++i)
279 outFile <<std::setprecision(20) << float_seed_table[i] << " ";
280 outFile << std::endl;
281 outFile << i_lag << " " << j_lag << std::endl;
282 outFile << std::setprecision(20) << carry << " " << count24 << std::endl;
283 outFile << luxury << " " << nskip << std::endl;
284 }
285 #endif
286 }
287
restoreStatus(const char filename[])288 void RanluxEngine::restoreStatus( const char filename[] )
289 {
290 std::ifstream inFile( filename, std::ios::in);
291 if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) {
292 std::cerr << " -- Engine state remains unchanged\n";
293 return;
294 }
295 if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) {
296 std::vector<unsigned long> v;
297 unsigned long xin;
298 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
299 inFile >> xin;
300 #ifdef TRACE_IO
301 std::cout << "ivec = " << ivec << " xin = " << xin << " ";
302 if (ivec%3 == 0) std::cout << "\n";
303 #endif
304 if (!inFile) {
305 inFile.clear(std::ios::badbit | inFile.rdstate());
306 std::cerr << "\nRanluxEngine state (vector) description improper."
307 << "\nrestoreStatus has failed."
308 << "\nInput stream is probably mispositioned now." << std::endl;
309 return;
310 }
311 v.push_back(xin);
312 }
313 getState(v);
314 return;
315 }
316
317 if (!inFile.bad() && !inFile.eof()) {
318 // inFile >> theSeed; removed -- encompased by possibleKeywordInput
319 for (int i=0; i<24; ++i)
320 inFile >> float_seed_table[i];
321 inFile >> i_lag; inFile >> j_lag;
322 inFile >> carry; inFile >> count24;
323 inFile >> luxury; inFile >> nskip;
324 }
325 }
326
showStatus() const327 void RanluxEngine::showStatus() const
328 {
329 std::cout << std::endl;
330 std::cout << "--------- Ranlux engine status ---------" << std::endl;
331 std::cout << " Initial seed = " << theSeed << std::endl;
332 std::cout << " float_seed_table[] = ";
333 for (int i=0; i<24; ++i)
334 std::cout << float_seed_table[i] << " ";
335 std::cout << std::endl;
336 std::cout << " i_lag = " << i_lag << ", j_lag = " << j_lag << std::endl;
337 std::cout << " carry = " << carry << ", count24 = " << count24 << std::endl;
338 std::cout << " luxury = " << luxury << " nskip = " << nskip << std::endl;
339 std::cout << "----------------------------------------" << std::endl;
340 }
341
flat()342 double RanluxEngine::flat() {
343
344 float next_random;
345 float uni;
346 int i;
347
348 uni = float_seed_table[j_lag] - float_seed_table[i_lag] - carry;
349 #ifdef TRACE_IO
350 if (flat_trace) {
351 std::cout << "float_seed_table[" << j_lag << "] = "
352 << float_seed_table[j_lag]
353 << " float_seed_table[" << i_lag << "] = " << float_seed_table[i_lag]
354 << " uni = " << uni << "\n";
355 std::cout << float_seed_table[j_lag]
356 << " - " << float_seed_table[i_lag]
357 << " - " << carry << " = "
358 << (double)float_seed_table[j_lag]
359 - (double) float_seed_table[i_lag] - (double)carry
360 << "\n";
361 }
362 #endif
363 if(uni < 0. ){
364 uni += 1.0;
365 carry = mantissa_bit_24();
366 }else{
367 carry = 0.;
368 }
369
370 float_seed_table[i_lag] = uni;
371 i_lag --;
372 j_lag --;
373 if(i_lag < 0) i_lag = 23;
374 if(j_lag < 0) j_lag = 23;
375
376 if( uni < mantissa_bit_12() ){
377 uni += mantissa_bit_24() * float_seed_table[j_lag];
378 if( uni == 0) uni = mantissa_bit_24() * mantissa_bit_24();
379 }
380 next_random = uni;
381 count24 ++;
382
383 // every 24th number generation, several random numbers are generated
384 // and wasted depending upon the luxury level.
385
386 if(count24 == 24 ){
387 count24 = 0;
388 #ifdef TRACE_IO
389 if (flat_trace) {
390 std::cout << "carry = " << carry << "\n";
391 }
392 #endif
393 for( i = 0; i != nskip ; i++){
394 uni = float_seed_table[j_lag] - float_seed_table[i_lag] - carry;
395 if(uni < 0. ){
396 uni += 1.0;
397 carry = mantissa_bit_24();
398 }else{
399 carry = 0.;
400 }
401 float_seed_table[i_lag] = uni;
402 #ifdef TRACE_IO
403 if (flat_trace) {
404 double xfst = float_seed_table[i_lag];
405 std::cout << "fst[" << i_lag << "] = "
406 << DoubConv::d2x(xfst) << "\n";
407 }
408 #endif
409 i_lag --;
410 j_lag --;
411 if(i_lag < 0)i_lag = 23;
412 if(j_lag < 0) j_lag = 23;
413 }
414 }
415 #ifdef TRACE_IO
416 if (flat_trace) {
417 std::cout << "next_random = " << next_random << "\n";
418 // flat_trace = false;
419 }
420 #endif
421 return (double) next_random;
422 }
423
flatArray(const int size,double * vect)424 void RanluxEngine::flatArray(const int size, double* vect)
425 {
426 float next_random;
427 float uni;
428 int i;
429 int index;
430
431 for (index=0; index<size; ++index) {
432 uni = float_seed_table[j_lag] - float_seed_table[i_lag] - carry;
433 if(uni < 0. ){
434 uni += 1.0;
435 carry = mantissa_bit_24();
436 }else{
437 carry = 0.;
438 }
439
440 float_seed_table[i_lag] = uni;
441 i_lag --;
442 j_lag --;
443 if(i_lag < 0) i_lag = 23;
444 if(j_lag < 0) j_lag = 23;
445
446 if( uni < mantissa_bit_12() ){
447 uni += mantissa_bit_24() * float_seed_table[j_lag];
448 if( uni == 0) uni = mantissa_bit_24() * mantissa_bit_24();
449 }
450 next_random = uni;
451 vect[index] = (double)next_random;
452 count24 ++;
453
454 // every 24th number generation, several random numbers are generated
455 // and wasted depending upon the luxury level.
456
457 if(count24 == 24 ){
458 count24 = 0;
459 for( i = 0; i != nskip ; i++){
460 uni = float_seed_table[j_lag] - float_seed_table[i_lag] - carry;
461 if(uni < 0. ){
462 uni += 1.0;
463 carry = mantissa_bit_24();
464 }else{
465 carry = 0.;
466 }
467 float_seed_table[i_lag] = uni;
468 i_lag --;
469 j_lag --;
470 if(i_lag < 0)i_lag = 23;
471 if(j_lag < 0) j_lag = 23;
472 }
473 }
474 }
475 }
476
operator double()477 RanluxEngine::operator double() {
478 return flat();
479 }
480
operator float()481 RanluxEngine::operator float() {
482 return float( flat() );
483 }
484
operator unsigned int()485 RanluxEngine::operator unsigned int() {
486 return ((unsigned int)(flat() * exponent_bit_32()) & 0xffffffff) |
487 (((unsigned int)(float_seed_table[i_lag]*exponent_bit_32())>>16) & 0xff);
488 // needed because Ranlux doesn't fill all bits of the double
489 // which therefore doesn't fill all bits of the integer.
490 }
491
put(std::ostream & os) const492 std::ostream & RanluxEngine::put ( std::ostream& os ) const
493 {
494 char beginMarker[] = "RanluxEngine-begin";
495 os << beginMarker << "\nUvec\n";
496 std::vector<unsigned long> v = put();
497 for (unsigned int i=0; i<v.size(); ++i) {
498 os << v[i] << "\n";
499 }
500 return os;
501 #ifdef REMOVED
502 char endMarker[] = "RanluxEngine-end";
503 int pr = os.precision(20);
504 os << " " << beginMarker << " ";
505 os << theSeed << "\n";
506 for (int i=0; i<24; ++i) {
507 os << float_seed_table[i] << "\n";
508 }
509 os << i_lag << " " << j_lag << "\n";
510 os << carry << " " << count24 << " ";
511 os << luxury << " " << nskip << "\n";
512 os << endMarker << "\n";
513 os.precision(pr);
514 return os;
515 #endif
516 }
517
put() const518 std::vector<unsigned long> RanluxEngine::put () const {
519 std::vector<unsigned long> v;
520 v.push_back (engineIDulong<RanluxEngine>());
521 #ifdef TRACE_IO
522 std::cout << "RanluxEngine put: ID is " << v[0] << "\n";
523 #endif
524 for (int i=0; i<24; ++i) {
525 v.push_back
526 (static_cast<unsigned long>(float_seed_table[i]/mantissa_bit_24()));
527 #ifdef TRACE_IO
528 std::cout << "v[" << i+1 << "] = " << v[i+1] <<
529 " float_seed_table[" << i << "] = " << float_seed_table[i] << "\n";
530 #endif
531 }
532 v.push_back(static_cast<unsigned long>(i_lag));
533 v.push_back(static_cast<unsigned long>(j_lag));
534 v.push_back(static_cast<unsigned long>(carry/mantissa_bit_24()));
535 v.push_back(static_cast<unsigned long>(count24));
536 v.push_back(static_cast<unsigned long>(luxury));
537 v.push_back(static_cast<unsigned long>(nskip));
538 #ifdef TRACE_IO
539 std::cout << "i_lag: " << v[25] << " j_lag: " << v[26]
540 << " carry: " << v[27] << "\n";
541 std::cout << "count24: " << v[28] << " luxury: " << v[29]
542 << " nskip: " << v[30] << "\n";
543 #endif
544 #ifdef TRACE_IO
545 flat_trace = true;
546 #endif
547 return v;
548 }
549
get(std::istream & is)550 std::istream & RanluxEngine::get ( std::istream& is )
551 {
552 char beginMarker [MarkerLen];
553 is >> std::ws;
554 is.width(MarkerLen); // causes the next read to the char* to be <=
555 // that many bytes, INCLUDING A TERMINATION \0
556 // (Stroustrup, section 21.3.2)
557 is >> beginMarker;
558 if (strcmp(beginMarker,"RanluxEngine-begin")) {
559 is.clear(std::ios::badbit | is.rdstate());
560 std::cerr << "\nInput stream mispositioned or"
561 << "\nRanluxEngine state description missing or"
562 << "\nwrong engine type found." << std::endl;
563 return is;
564 }
565 return getState(is);
566 }
567
beginTag()568 std::string RanluxEngine::beginTag ( ) {
569 return "RanluxEngine-begin";
570 }
571
getState(std::istream & is)572 std::istream & RanluxEngine::getState ( std::istream& is )
573 {
574 if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) {
575 std::vector<unsigned long> v;
576 unsigned long uu;
577 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
578 is >> uu;
579 if (!is) {
580 is.clear(std::ios::badbit | is.rdstate());
581 std::cerr << "\nRanluxEngine state (vector) description improper."
582 << "\ngetState() has failed."
583 << "\nInput stream is probably mispositioned now." << std::endl;
584 return is;
585 }
586 v.push_back(uu);
587 #ifdef TRACE_IO
588 std::cout << "RanluxEngine::getState -- v[" << v.size()-1
589 << "] = " << v[v.size()-1] << "\n";
590 #endif
591 }
592 getState(v);
593 return (is);
594 }
595
596 // is >> theSeed; Removed, encompassed by possibleKeywordInput()
597
598 char endMarker [MarkerLen];
599 for (int i=0; i<24; ++i) {
600 is >> float_seed_table[i];
601 }
602 is >> i_lag; is >> j_lag;
603 is >> carry; is >> count24;
604 is >> luxury; is >> nskip;
605 is >> std::ws;
606 is.width(MarkerLen);
607 is >> endMarker;
608 if (strcmp(endMarker,"RanluxEngine-end")) {
609 is.clear(std::ios::badbit | is.rdstate());
610 std::cerr << "\nRanluxEngine state description incomplete."
611 << "\nInput stream is probably mispositioned now." << std::endl;
612 return is;
613 }
614 return is;
615 }
616
get(const std::vector<unsigned long> & v)617 bool RanluxEngine::get (const std::vector<unsigned long> & v) {
618 if ((v[0] & 0xffffffffUL) != engineIDulong<RanluxEngine>()) {
619 std::cerr <<
620 "\nRanluxEngine get:state vector has wrong ID word - state unchanged\n";
621 return false;
622 }
623 return getState(v);
624 }
625
getState(const std::vector<unsigned long> & v)626 bool RanluxEngine::getState (const std::vector<unsigned long> & v) {
627 if (v.size() != VECTOR_STATE_SIZE ) {
628 std::cerr <<
629 "\nRanluxEngine get:state vector has wrong length - state unchanged\n";
630 return false;
631 }
632 for (int i=0; i<24; ++i) {
633 float_seed_table[i] = v[i+1]*mantissa_bit_24();
634 #ifdef TRACE_IO
635 std::cout <<
636 "float_seed_table[" << i << "] = " << float_seed_table[i] << "\n";
637 #endif
638 }
639 i_lag = v[25];
640 j_lag = v[26];
641 carry = v[27]*mantissa_bit_24();
642 count24 = v[28];
643 luxury = v[29];
644 nskip = v[30];
645 #ifdef TRACE_IO
646 std::cout << "i_lag: " << i_lag << " j_lag: " << j_lag
647 << " carry: " << carry << "\n";
648 std::cout << "count24: " << count24 << " luxury: " << luxury
649 << " nskip: " << nskip << "\n";
650
651 #endif
652 #ifdef TRACE_IO
653 flat_trace = true;
654 #endif
655 return true;
656 }
657
658 } // namespace CLHEP
659