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

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