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

RandEngine.cc
Go to the documentation of this file.
1 // $Id: RandEngine.cc,v 1.8 2010/06/16 17:24:53 garren Exp $
2 // -*- C++ -*-
3 //
4 // -----------------------------------------------------------------------
5 // HEP Random
6 // --- RandEngine ---
7 // class implementation file
8 // -----------------------------------------------------------------------
9 // This file is part of Geant4 (simulation toolkit for HEP).
10 
11 // =======================================================================
12 // Gabriele Cosmo - Created: 5th September 1995
13 // - Minor corrections: 31st October 1996
14 // - Added methods for engine status: 19th November 1996
15 // - mx is initialised to RAND_MAX: 2nd April 1997
16 // - Fixed bug in setSeeds(): 15th September 1997
17 // - Private copy constructor and operator=: 26th Feb 1998
18 // J.Marraffino - Added stream operators and related constructor.
19 // Added automatic seed selection from seed table and
20 // engine counter. Removed RAND_MAX and replaced by
21 // std::pow(0.5,32.). Flat() returns now 32 bits values
22 // obtained by concatenation: 15th Feb 1998
23 // Ken Smith - Added conversion operators: 6th Aug 1998
24 // J. Marraffino - Remove dependence on hepString class 13 May 1999
25 // M. Fischler - Rapaired bug that in flat() that relied on rand() to
26 // deliver 15-bit results. This bug was causing flat()
27 // on Linux systems to yield randoms with mean of 5/8(!)
28 // - Modified algorithm such that on systems with 31-bit rand()
29 // it will call rand() only once instead of twice. Feb 2004
30 // M. Fischler - Modified the general-case template for RandEngineBuilder
31 // such that when RAND_MAX is an unexpected value the routine
32 // will still deliver a sensible flat() random.
33 // M. Fischler - Methods for distrib. instance save/restore 12/8/04
34 // M. Fischler - split get() into tag validation and
35 // getState() for anonymous restores 12/27/04
36 // M. Fischler - put/get for vectors of ulongs 3/14/05
37 // M. Fischler - State-saving using only ints, for portability 4/12/05
38 //
39 // =======================================================================
40 
41 #include "CLHEP/Random/defs.h"
42 #include "CLHEP/Random/RandEngine.h"
43 #include "CLHEP/Random/Random.h"
44 #include "CLHEP/Random/engineIDulong.h"
45 #include <string.h> // for strcmp
46 #include <cstdlib> // for int()
47 
48 namespace CLHEP {
49 
50 #ifdef NOTDEF
51 // The way to test for proper behavior of the RandEngineBuilder
52 // for arbitrary RAND_MAX, on a system where RAND_MAX is some
53 // fixed specialized value and rand() behaves accordingly, is
54 // to set up a fake RAND_MAX and a fake version of rand()
55 // by enabling this block.
56 #undef RAND_MAX
57 #define RAND_MAX 9382956
58 #include "CLHEP/Random/MTwistEngine.h"
59 #include "CLHEP/Random/RandFlat.h"
60 MTwistEngine * fakeFlat = new MTwistEngine;
61 RandFlat rflat (fakeFlat, 0, RAND_MAX+1);
62 int rand() { return (int)rflat(); }
63 #endif
64 
65 static const int MarkerLen = 64; // Enough room to hold a begin or end marker.
66 
67 // number of instances with automatic seed selection
68 int RandEngine::numEngines = 0;
69 
70 // Maximum index into the seed table
71 int RandEngine::maxIndex = 215;
72 
73 std::string RandEngine::name() const {return "RandEngine";}
74 
77 {
78  setSeed(seed,0);
79  setSeeds(&theSeed,0);
80  seq = 0;
81 }
82 
83 RandEngine::RandEngine(int rowIndex, int colIndex)
85 {
86  long seeds[2];
87  long seed;
88 
89  int cycle = std::abs(int(rowIndex/maxIndex));
90  int row = std::abs(int(rowIndex%maxIndex));
91  int col = std::abs(int(colIndex%2));
92  long mask = ((cycle & 0x000007ff) << 20 );
93  HepRandom::getTheTableSeeds( seeds, row );
94  seed = (seeds[col])^mask;
95  setSeed(seed,0);
96  setSeeds(&theSeed,0);
97  seq = 0;
98 }
99 
101 : HepRandomEngine()
102 {
103  long seeds[2];
104  long seed;
105  int cycle,curIndex;
106 
107  cycle = std::abs(int(numEngines/maxIndex));
108  curIndex = std::abs(int(numEngines%maxIndex));
109  numEngines += 1;
110  long mask = ((cycle & 0x007fffff) << 8);
111  HepRandom::getTheTableSeeds( seeds, curIndex );
112  seed = seeds[0]^mask;
113  setSeed(seed,0);
114  setSeeds(&theSeed,0);
115  seq = 0;
116 }
117 
118 RandEngine::RandEngine(std::istream& is)
119 : HepRandomEngine()
120 {
121  is >> *this;
122 }
123 
125 
126 void RandEngine::setSeed(long seed, int)
127 {
128  theSeed = seed;
129  srand( int(seed) );
130  seq = 0;
131 }
132 
133 void RandEngine::setSeeds(const long* seeds, int)
134 {
135  setSeed(seeds ? *seeds : 19780503L, 0);
136  theSeeds = seeds;
137 }
138 
139 void RandEngine::saveStatus( const char filename[] ) const
140 {
141  std::ofstream outFile( filename, std::ios::out ) ;
142 
143  if (!outFile.bad()) {
144  outFile << "Uvec\n";
145  std::vector<unsigned long> v = put();
146  #ifdef TRACE_IO
147  std::cout << "Result of v = put() is:\n";
148  #endif
149  for (unsigned int i=0; i<v.size(); ++i) {
150  outFile << v[i] << "\n";
151  #ifdef TRACE_IO
152  std::cout << v[i] << " ";
153  if (i%6==0) std::cout << "\n";
154  #endif
155  }
156  #ifdef TRACE_IO
157  std::cout << "\n";
158  #endif
159  }
160 #ifdef REMOVED
161  if (!outFile.bad()) {
162  outFile << theSeed << std::endl;
163  outFile << seq << std::endl;
164  }
165 #endif
166 }
167 
168 void RandEngine::restoreStatus( const char filename[] )
169 {
170  // The only way to restore the status of RandEngine is to
171  // keep track of the number of shooted random sequences, reset
172  // the engine and re-shoot them again. The Rand algorithm does
173  // not provide any way of getting its internal status.
174 
175  std::ifstream inFile( filename, std::ios::in);
176  if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) {
177  std::cout << " -- Engine state remains unchanged\n";
178  return;
179  }
180  if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) {
181  std::vector<unsigned long> v;
182  unsigned long xin;
183  for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
184  inFile >> xin;
185  #ifdef TRACE_IO
186  std::cout << "ivec = " << ivec << " xin = " << xin << " ";
187  if (ivec%3 == 0) std::cout << "\n";
188  #endif
189  if (!inFile) {
190  inFile.clear(std::ios::badbit | inFile.rdstate());
191  std::cerr << "\nRandEngine state (vector) description improper."
192  << "\nrestoreStatus has failed."
193  << "\nInput stream is probably mispositioned now." << std::endl;
194  return;
195  }
196  v.push_back(xin);
197  }
198  getState(v);
199  return;
200  }
201 
202  long count;
203 
204  if (!inFile.bad() && !inFile.eof()) {
205 // inFile >> theSeed; removed -- encompased by possibleKeywordInput
206  inFile >> count;
207  setSeed(theSeed,0);
208  seq = 0;
209  while (seq < count) flat();
210  }
211 }
212 
214 {
215  std::cout << std::endl;
216  std::cout << "---------- Rand engine status ----------" << std::endl;
217  std::cout << " Initial seed = " << theSeed << std::endl;
218  std::cout << " Shooted sequences = " << seq << std::endl;
219  std::cout << "----------------------------------------" << std::endl;
220 }
221 
222 // ====================================================
223 // Implementation of flat() (along with needed helpers)
224 // ====================================================
225 
226 // Here we set up such that **at compile time**, the compiler decides based on
227 // RAND_MAX how to form a random double with 32 leading random bits, using
228 // one or two calls to rand(). Some of the lowest order bits of 32 are allowed
229 // to be as weak as mere XORs of some higher bits, but not to be always fixed.
230 //
231 // The method decision is made at compile time, rather than using a run-time
232 // if on the value of RAND_MAX. Although it is possible to cope with arbitrary
233 // values of RAND_MAX of the form 2**N-1, with the same efficiency, the
234 // template techniques needed would look mysterious and perhaps over-stress
235 // some older compilers. We therefore only treat RAND_MAX = 2**15-1 (as on
236 // most older systems) and 2**31-1 (as on the later Linux systems) as special
237 // and super-efficient cases. We detect any different value, and use an
238 // algorithm which is correct even if RAND_MAX is not one less than a power
239 // of 2.
240 
241  template <int> struct RandEngineBuilder { // RAND_MAX any arbitrary value
242  static unsigned int thirtyTwoRandomBits(long& seq) {
243 
244  static bool prepared = false;
245  static unsigned int iT;
246  static unsigned int iK;
247  static unsigned int iS;
248  static int iN;
249  static double fS;
250  static double fT;
251 
252  if ( (RAND_MAX >> 31) > 0 )
253  {
254  // Here, we know that integer arithmetic is 64 bits.
255  if ( !prepared ) {
256  iS = (unsigned long)RAND_MAX + 1;
257  iK = 1;
258 // int StoK = S;
259  int StoK = iS;
260  // The two statements below are equivalent, but some compilers
261  // are getting too smart and complain about the first statement.
262  //if ( (RAND_MAX >> 32) == 0) {
263  if( (unsigned long) (RAND_MAX) <= (( (1uL) << 31 ) - 1 ) ) {
264  iK = 2;
265 // StoK = S*S;
266  StoK = iS*iS;
267  }
268  int m;
269  for ( m = 0; m < 64; ++m ) {
270  StoK >>= 1;
271  if (StoK == 0) break;
272  }
273  iT = 1 << m;
274  prepared = true;
275  }
276  int v = 0;
277  do {
278  for ( int i = 0; i < iK; ++i ) {
279  v = iS*v+rand(); ++seq;
280  }
281  } while (v < iT);
282  return v & 0xFFFFFFFF;
283 
284  }
285 
286  else if ( (RAND_MAX >> 26) == 0 )
287  {
288  // Here, we know it is safe to work in terms of doubles without loss
289  // of precision, but we have no idea how many randoms we will need to
290  // generate 32 bits.
291  if ( !prepared ) {
292  fS = (unsigned long)RAND_MAX + 1;
293  double twoTo32 = std::ldexp(1.0,32);
294  double StoK = fS;
295  for ( iK = 1; StoK < twoTo32; StoK *= fS, iK++ ) { }
296  int m;
297  fT = 1.0;
298  for ( m = 0; m < 64; ++m ) {
299  StoK *= .5;
300  if (StoK < 1.0) break;
301  fT *= 2.0;
302  }
303  prepared = true;
304  }
305  double v = 0;
306  do {
307  for ( int i = 0; i < iK; ++i ) {
308  v = fS*v+rand(); ++seq;
309  }
310  } while (v < fT);
311  return ((unsigned int)v) & 0xFFFFFFFF;
312 
313  }
314  else
315  {
316  // Here, we know that 16 random bits are available from each of
317  // two random numbers.
318  if ( !prepared ) {
319  iS = (unsigned long)RAND_MAX + 1;
320  int SshiftN = iS;
321  for (iN = 0; SshiftN > 1; SshiftN >>= 1, iN++) { }
322  iN -= 17;
323  prepared = true;
324  }
325  unsigned int x1, x2;
326  do { x1 = rand(); ++seq;} while (x1 < (1<<16) );
327  do { x2 = rand(); ++seq;} while (x2 < (1<<16) );
328  return x1 | (x2 << 16);
329  }
330 
331  }
332 };
333 
334 template <> struct RandEngineBuilder<2147483647> { // RAND_MAX = 2**31 - 1
335  inline static unsigned int thirtyTwoRandomBits(long& seq) {
336  unsigned int x = rand() << 1; ++seq; // bits 31-1
337  x ^= ( (x>>23) ^ (x>>7) ) ^1; // bit 0 (weakly pseudo-random)
338  return x & 0xFFFFFFFF; // mask in case int is 64 bits
339  }
340 };
341 
342 
343 template <> struct RandEngineBuilder<32767> { // RAND_MAX = 2**15 - 1
344  inline static unsigned int thirtyTwoRandomBits(long& seq) {
345  unsigned int x = rand() << 17; ++seq; // bits 31-17
346  x ^= rand() << 2; ++seq; // bits 16-2
347  x ^= ( (x>>23) ^ (x>>7) ) ^3; // bits 1-0 (weakly pseudo-random)
348  return x & 0xFFFFFFFF; // mask in case int is 64 bits
349  }
350 };
351 
353 {
354  double r;
356  } while ( r == 0 );
357  return r/4294967296.0;
358 }
359 
360 void RandEngine::flatArray(const int size, double* vect)
361 {
362  int i;
363 
364  for (i=0; i<size; ++i)
365  vect[i]=flat();
366 }
367 
368 RandEngine::operator unsigned int() {
370 }
371 
372 std::ostream & RandEngine::put ( std::ostream& os ) const
373 {
374  char beginMarker[] = "RandEngine-begin";
375  char endMarker[] = "RandEngine-end";
376 
377  os << " " << beginMarker << "\n";
378  os << theSeed << " " << seq << " ";
379  os << endMarker << "\n";
380  return os;
381 }
382 
383 std::vector<unsigned long> RandEngine::put () const {
384  std::vector<unsigned long> v;
385  v.push_back (engineIDulong<RandEngine>());
386  v.push_back(static_cast<unsigned long>(theSeed));
387  v.push_back(static_cast<unsigned long>(seq));
388  return v;
389 }
390 
391 std::istream & RandEngine::get ( std::istream& is )
392 {
393  // The only way to restore the status of RandEngine is to
394  // keep track of the number of shooted random sequences, reset
395  // the engine and re-shoot them again. The Rand algorithm does
396  // not provide any way of getting its internal status.
397  char beginMarker [MarkerLen];
398  is >> std::ws;
399  is.width(MarkerLen); // causes the next read to the char* to be <=
400  // that many bytes, INCLUDING A TERMINATION \0
401  // (Stroustrup, section 21.3.2)
402  is >> beginMarker;
403  if (strcmp(beginMarker,"RandEngine-begin")) {
404  is.clear(std::ios::badbit | is.rdstate());
405  std::cout << "\nInput stream mispositioned or"
406  << "\nRandEngine state description missing or"
407  << "\nwrong engine type found." << std::endl;
408  return is;
409  }
410  return getState(is);
411 }
412 
413 std::string RandEngine::beginTag ( ) {
414  return "RandEngine-begin";
415 }
416 
417 std::istream & RandEngine::getState ( std::istream& is )
418 {
419  if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) {
420  std::vector<unsigned long> v;
421  unsigned long uu;
422  for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
423  is >> uu;
424  if (!is) {
425  is.clear(std::ios::badbit | is.rdstate());
426  std::cerr << "\nRandEngine state (vector) description improper."
427  << "\ngetState() has failed."
428  << "\nInput stream is probably mispositioned now." << std::endl;
429  return is;
430  }
431  v.push_back(uu);
432  }
433  getState(v);
434  return (is);
435  }
436 
437 // is >> theSeed; Removed, encompassed by possibleKeywordInput()
438 
439  char endMarker [MarkerLen];
440  long count;
441  is >> count;
442  is >> std::ws;
443  is.width(MarkerLen);
444  is >> endMarker;
445  if (strcmp(endMarker,"RandEngine-end")) {
446  is.clear(std::ios::badbit | is.rdstate());
447  std::cerr << "\nRandEngine state description incomplete."
448  << "\nInput stream is probably mispositioned now." << std::endl;
449  return is;
450  }
451  setSeed(theSeed,0);
452  while (seq < count) flat();
453  return is;
454 }
455 
456 bool RandEngine::get (const std::vector<unsigned long> & v) {
457  if ((v[0] & 0xffffffffUL) != engineIDulong<RandEngine>()) {
458  std::cerr <<
459  "\nRandEngine get:state vector has wrong ID word - state unchanged\n";
460  return false;
461  }
462  return getState(v);
463 }
464 
465 bool RandEngine::getState (const std::vector<unsigned long> & v) {
466  if (v.size() != VECTOR_STATE_SIZE ) {
467  std::cerr <<
468  "\nRandEngine get:state vector has wrong length - state unchanged\n";
469  return false;
470  }
471  theSeed = v[1];
472  int count = v[2];
473  setSeed(theSeed,0);
474  while (seq < count) flat();
475  return true;
476 }
477 } // namespace CLHEP
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
static std::string engineName()
virtual std::istream & getState(std::istream &is)
Definition: RandEngine.cc:417
void flatArray(const int size, double *vect)
Definition: RandEngine.cc:360
static unsigned int thirtyTwoRandomBits(long &seq)
Definition: RandEngine.cc:344
std::string name() const
Definition: RandEngine.cc:73
void setSeeds(const long *seeds, int dum=0)
Definition: RandEngine.cc:133
void saveStatus(const char filename[]="Rand.conf") const
Definition: RandEngine.cc:139
std::vector< unsigned long > put() const
Definition: RandEngine.cc:383
static const unsigned int VECTOR_STATE_SIZE
void restoreStatus(const char filename[]="Rand.conf")
Definition: RandEngine.cc:168
static unsigned int thirtyTwoRandomBits(long &seq)
Definition: RandEngine.cc:242
void setSeed(long seed, int dum=0)
Definition: RandEngine.cc:126
virtual ~RandEngine()
Definition: RandEngine.cc:124
static bool checkFile(std::istream &file, const std::string &filename, const std::string &classname, const std::string &methodname)
Definition: RandomEngine.cc:46
static void getTheTableSeeds(long *seeds, int index)
Definition: Random.cc:152
static std::string beginTag()
Definition: RandEngine.cc:413
virtual std::istream & get(std::istream &is)
Definition: RandEngine.cc:391
void showStatus() const
Definition: RandEngine.cc:213
static unsigned int thirtyTwoRandomBits(long &seq)
Definition: RandEngine.cc:335