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