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