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 
52 namespace CLHEP {
53 
54 
55 static const int MarkerLen = 64; // Enough room to hold a begin or end marker.
56 
57 std::string RanluxEngine::name() const {return "RanluxEngine";}
58 
59 // Number of instances with automatic seed selection
60 int RanluxEngine::numEngines = 0;
61 
62 // Maximum index into the seed table
63 int RanluxEngine::maxIndex = 215;
64 
65 RanluxEngine::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 
100 RanluxEngine::RanluxEngine(int rowIndex, int colIndex, int lux)
101 : HepRandomEngine()
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 
121 RanluxEngine::RanluxEngine( std::istream& is )
122 : HepRandomEngine()
123 {
124  is >> *this;
125 }
126 
128 
129 void 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 
185 void 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 
251 void 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 
284 void 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 
420 void 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 
473 RanluxEngine::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 
480 std::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 
506 std::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 
538 std::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 
556 std::string RanluxEngine::beginTag ( ) {
557  return "RanluxEngine-begin";
558 }
559 
560 std::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 
605 bool 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 
614 bool 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
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
#define double(obj)
Definition: excDblThrow.cc:32
void setSeed(long seed, int lux=3)
void showStatus() const
virtual std::istream & getState(std::istream &is)
void flatArray(const int size, double *vect)
static double mantissa_bit_12()
static double mantissa_bit_24()
void setSeeds(const long *seeds, int lux=3)
virtual std::istream & get(std::istream &is)
std::string name() const
Definition: RanluxEngine.cc:57
static bool checkFile(std::istream &file, const std::string &filename, const std::string &classname, const std::string &methodname)
Definition: RandomEngine.cc:46
static std::string d2x(double d)
Definition: DoubConv.cc:78
static void getTheTableSeeds(long *seeds, int index)
Definition: Random.cc:152
std::vector< unsigned long > put() const
static const unsigned int VECTOR_STATE_SIZE
void saveStatus(const char filename[]="Ranlux.conf") const
static std::string beginTag()
void restoreStatus(const char filename[]="Ranlux.conf")