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 );
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 
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
CLHEP::RandEngine::restoreStatus
void restoreStatus(const char filename[]="Rand.conf")
Definition: RandEngine.cc:168
CLHEP::HepRandomEngine
Definition: Matrix/CLHEP/Random/RandomEngine.h:55
CLHEP::RandEngineBuilder::thirtyTwoRandomBits
static unsigned int thirtyTwoRandomBits(long &seq)
Definition: RandEngine.cc:242
CLHEP::HepRandomEngine::theSeed
long theSeed
Definition: Matrix/CLHEP/Random/RandomEngine.h:144
CLHEP::RandEngineBuilder< 2147483647 >::thirtyTwoRandomBits
static unsigned int thirtyTwoRandomBits(long &seq)
Definition: RandEngine.cc:335
CLHEP::RandEngine::~RandEngine
virtual ~RandEngine()
Definition: RandEngine.cc:124
CLHEP::HepRandom::getTheTableSeeds
static void getTheTableSeeds(long *seeds, int index)
Definition: Random.cc:152
CLHEP::RandEngineBuilder< 32767 >::thirtyTwoRandomBits
static unsigned int thirtyTwoRandomBits(long &seq)
Definition: RandEngine.cc:344
is
HepRotation and so forth isNear() norm2() rectify() static Rotation row1 row4(To avoid bloat in the code pulled in for programs which don 't use all these features, we split the implementation .cc files. Only isNear() goes into the original Rotation.cc) --------------------------------------- HepAxisAngle and HepEulerAngles classes --------------------------------------- These classes are very useful and simple structures for holding the result of a nice intuituve decomposition of a rotation there is no longer much content in the distinct ZOOM PhysicsVectors library The only content left in the library is the object files representing the various Exception objects When we build the CLHEP classes for the ZOOM we will set up so as to use ZOOM SpaceVector is(but we can disable namespace usage and most of our users do so at this point). What I do is leave Hep3Vector in the global namespace
CLHEP::HepRandomEngine::theSeeds
const long * theSeeds
Definition: Matrix/CLHEP/Random/RandomEngine.h:145
CLHEP::RandEngine::beginTag
static std::string beginTag()
Definition: RandEngine.cc:413
CLHEP::RandEngine::saveStatus
void saveStatus(const char filename[]="Rand.conf") const
Definition: RandEngine.cc:139
CLHEP::RandEngine::setSeeds
void setSeeds(const long *seeds, int dum=0)
Definition: RandEngine.cc:133
size
user code seldom needs to call this function directly ZMerrno whether or not they are still recorded ZMerrno size() Return the(integer) number of ZMthrow 'n exceptions currently recorded. 5) ZMerrno.clear() Set an internal counter to zero. This counter is available(see next function) to user code to track ZMthrow 'n exceptions that have occurred during any arbitrary time interval. 6) ZMerrno.countSinceCleared() Return the(integer) number of ZMthrow 'n exceptions that have been recorded via ZMerrno.write()
CLHEP
Definition: ClhepVersion.h:13
v
they are gone ZOOM Features Discontinued The following features of the ZOOM package were felt to be extreme overkill These have been after checking that no existing user code was utilizing as in SpaceVector v
Definition: keyMergeIssues.doc:324
CLHEP::RandEngine::name
std::string name() const
Definition: RandEngine.cc:73
CLHEP::RandEngine::RandEngine
RandEngine()
Definition: RandEngine.cc:100
CLHEP::RandEngine::put
std::vector< unsigned long > put() const
Definition: RandEngine.cc:383
CLHEP::RandEngine::get
virtual std::istream & get(std::istream &is)
Definition: RandEngine.cc:391
CLHEP::RandEngine::engineName
static std::string engineName()
Definition: Matrix/CLHEP/Random/RandEngine.h:98
CLHEP::RandEngine::flatArray
void flatArray(const int size, double *vect)
Definition: RandEngine.cc:360
CLHEP::RandEngineBuilder
Definition: RandEngine.cc:241
CLHEP::possibleKeywordInput
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
Definition: Matrix/CLHEP/Random/RandomEngine.h:168
CLHEP::RandEngine::setSeed
void setSeed(long seed, int dum=0)
Definition: RandEngine.cc:126
seeds
Technical Maintenance Note for CLHEP Random Consequences of seeding JamesRandom with positive seed values greater than In the source code JamesRandom The usual way of seeding a generator is via the default which makes use of the table of seeds(with some trickery to ensure that the values won 't repeat after the table rows are exhausted). The trickery preserves the fact that sees are never negative(because the table values are never negative
i
long i
Definition: JamesRandomSeeding.txt:27
CLHEP::RandEngine::getState
virtual std::istream & getState(std::istream &is)
Definition: RandEngine.cc:417
CLHEP::RandEngine::VECTOR_STATE_SIZE
static const unsigned int VECTOR_STATE_SIZE
Definition: Matrix/CLHEP/Random/RandEngine.h:104
CLHEP::RandEngine::showStatus
void showStatus() const
Definition: RandEngine.cc:213
CLHEP::HepRandomEngine::checkFile
static bool checkFile(std::istream &file, const std::string &filename, const std::string &classname, const std::string &methodname)
Definition: RandomEngine.cc:46
in
it has advantages For I leave the ZMthrows in
Definition: keyMergeIssues.doc:62
count
user code seldom needs to call this function directly ZMerrno count() Return the(integer) number of ZMthrow 'n exceptions ever recorded via ZMerrno.write()
x
any side effects of that construction would occur twice The semantics of throw x
Definition: whyZMthrowRethrows.txt:37
CLHEP::RandEngine::flat
double flat()
Definition: RandEngine.cc:352