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