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

RanecuEngine.cc
Go to the documentation of this file.
1 // $Id: RanecuEngine.cc,v 1.7 2010/07/20 18:06:02 garren Exp $
2 // -*- C++ -*-
3 //
4 // -----------------------------------------------------------------------
5 // HEP Random
6 // --- RanecuEngine ---
7 // class implementation file
8 // -----------------------------------------------------------------------
9 // This file is part of Geant4 (simulation toolkit for HEP).
10 //
11 // RANECU Random Engine - algorithm originally written in FORTRAN77
12 // as part of the MATHLIB HEP library.
13 
14 // =======================================================================
15 // Gabriele Cosmo - Created - 2nd February 1996
16 // - Minor corrections: 31st October 1996
17 // - Added methods for engine status: 19th November 1996
18 // - Added abs for setting seed index: 11th July 1997
19 // - Modified setSeeds() to handle default index: 16th Oct 1997
20 // - setSeed() now resets the engine status to the original
21 // values in the static table of HepRandom: 19th Mar 1998
22 // J.Marraffino - Added stream operators and related constructor.
23 // Added automatic seed selection from seed table and
24 // engine counter: 16th Feb 1998
25 // Ken Smith - Added conversion operators: 6th Aug 1998
26 // J. Marraffino - Remove dependence on hepString class 13 May 1999
27 // M. Fischler - Add endl to the end of saveStatus 10 Apr 2001
28 // M. Fischler - In restore, checkFile for file not found 03 Dec 2004
29 // M. Fischler - Methods for distrib. instance save/restore 12/8/04
30 // M. Fischler - split get() into tag validation and
31 // getState() for anonymous restores 12/27/04
32 // M. Fischler - put/get for vectors of ulongs 3/14/05
33 // M. Fischler - State-saving using only ints, for portability 4/12/05
34 // M. Fischler - Modify ctor and setSeed to utilize all info provided
35 // and avoid coincidence of same state from different
36 // seeds 6/22/10
37 //
38 // =======================================================================
39 
40 #include "CLHEP/Random/defs.h"
41 #include "CLHEP/Random/Random.h"
42 #include "CLHEP/Random/RanecuEngine.h"
43 #include "CLHEP/Random/engineIDulong.h"
44 #include <string.h> // for strcmp
45 #include <cmath>
46 #include <cstdlib>
47 
48 namespace CLHEP {
49 
50 static const int MarkerLen = 64; // Enough room to hold a begin or end marker.
51 
52 static const double prec = 4.6566128E-10;
53 
54 std::string RanecuEngine::name() const {return "RanecuEngine";}
55 
56 void RanecuEngine::further_randomize (int seq1, int col, int index, int modulus)
57 {
58  table[seq1][col] -= (index&0x3FFFFFFF);
59  while (table[seq1][col] <= 0) table[seq1][col] += (modulus-1);
60 } // mf 6/22/10
61 
62 // Number of instances with automatic seed selection
63 int RanecuEngine::numEngines = 0;
64 
67 {
68  int cycle = std::abs(int(numEngines/maxSeq));
69  seq = std::abs(int(numEngines%maxSeq));
70  numEngines += 1;
71  theSeed = seq;
72  long mask = ((cycle & 0x007fffff) << 8);
73  for (int i=0; i<2; ++i) {
74  for (int j=0; j<maxSeq; ++j) {
76  table[j][i] ^= mask;
77  }
78  }
79  theSeeds = &table[seq][0];
80 }
81 
84 {
85  int cycle = std::abs(int(index/maxSeq));
86  seq = std::abs(int(index%maxSeq));
87  theSeed = seq;
88  long mask = ((cycle & 0x000007ff) << 20);
89  for (int j=0; j<maxSeq; ++j) {
91  table[j][0] ^= mask;
92  table[j][1] ^= mask;
93  }
94  theSeeds = &table[seq][0];
95  further_randomize (seq, 0, index, shift1); // mf 6/22/10
96 }
97 
100 {
101  is >> *this;
102 }
103 
105 
106 void RanecuEngine::setSeed(long index, int dum)
107 {
108  seq = std::abs(int(index%maxSeq));
109  theSeed = seq;
110  HepRandom::getTheTableSeeds(table[seq],seq);
111  theSeeds = &table[seq][0];
112  further_randomize (seq, 0, index, shift1); // mf 6/22/10
113  further_randomize (seq, 1, dum, shift2); // mf 6/22/10
114 }
115 
116 void RanecuEngine::setSeeds(const long* seeds, int pos)
117 {
118  if (pos != -1) {
119  seq = std::abs(int(pos%maxSeq));
120  theSeed = seq;
121  }
122  // only positive seeds are allowed
123  table[seq][0] = std::abs(seeds[0])%shift1;
124  table[seq][1] = std::abs(seeds[1])%shift2;
125  theSeeds = &table[seq][0];
126 }
127 
128 void RanecuEngine::setIndex(long index)
129 {
130  seq = std::abs(int(index%maxSeq));
131  theSeed = seq;
132  theSeeds = &table[seq][0];
133 }
134 
135 void RanecuEngine::saveStatus( const char filename[] ) const
136 {
137  std::ofstream outFile( filename, std::ios::out ) ;
138 
139  if (!outFile.bad()) {
140  outFile << "Uvec\n";
141  std::vector<unsigned long> v = put();
142  #ifdef TRACE_IO
143  std::cout << "Result of v = put() is:\n";
144  #endif
145  for (unsigned int i=0; i<v.size(); ++i) {
146  outFile << v[i] << "\n";
147  #ifdef TRACE_IO
148  std::cout << v[i] << " ";
149  if (i%6==0) std::cout << "\n";
150  #endif
151  }
152  #ifdef TRACE_IO
153  std::cout << "\n";
154  #endif
155  }
156 #ifdef REMOVED
157  if (!outFile.bad()) {
158  outFile << theSeed << std::endl;
159  for (int i=0; i<2; ++i)
160  outFile << table[theSeed][i] << " ";
161  outFile << std::endl;
162  }
163 #endif
164 }
165 
166 void RanecuEngine::restoreStatus( const char filename[] )
167 {
168  std::ifstream inFile( filename, std::ios::in);
169  if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) {
170  std::cerr << " -- Engine state remains unchanged\n";
171  return;
172  }
173  if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) {
174  std::vector<unsigned long> v;
175  unsigned long xin;
176  for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
177  inFile >> xin;
178  #ifdef TRACE_IO
179  std::cout << "ivec = " << ivec << " xin = " << xin << " ";
180  if (ivec%3 == 0) std::cout << "\n";
181  #endif
182  if (!inFile) {
183  inFile.clear(std::ios::badbit | inFile.rdstate());
184  std::cerr << "\nJamesRandom state (vector) description improper."
185  << "\nrestoreStatus has failed."
186  << "\nInput stream is probably mispositioned now." << std::endl;
187  return;
188  }
189  v.push_back(xin);
190  }
191  getState(v);
192  return;
193  }
194 
195  if (!inFile.bad() && !inFile.eof()) {
196 // inFile >> theSeed; removed -- encompased by possibleKeywordInput
197  for (int i=0; i<2; ++i)
198  inFile >> table[theSeed][i];
199  seq = int(theSeed);
200  }
201 }
202 
204 {
205  std::cout << std::endl;
206  std::cout << "--------- Ranecu engine status ---------" << std::endl;
207  std::cout << " Initial seed (index) = " << theSeed << std::endl;
208  std::cout << " Current couple of seeds = "
209  << table[theSeed][0] << ", "
210  << table[theSeed][1] << std::endl;
211  std::cout << "----------------------------------------" << std::endl;
212 }
213 
215 {
216  const int index = seq;
217  long seed1 = table[index][0];
218  long seed2 = table[index][1];
219 
220  int k1 = (int)(seed1/ecuyer_b);
221  int k2 = (int)(seed2/ecuyer_e);
222 
223  seed1 = ecuyer_a*(seed1-k1*ecuyer_b)-k1*ecuyer_c;
224  if (seed1 < 0) seed1 += shift1;
225  seed2 = ecuyer_d*(seed2-k2*ecuyer_e)-k2*ecuyer_f;
226  if (seed2 < 0) seed2 += shift2;
227 
228  table[index][0] = seed1;
229  table[index][1] = seed2;
230 
231  long diff = seed1-seed2;
232 
233  if (diff <= 0) diff += (shift1-1);
234  return (double)(diff*prec);
235 }
236 
237 void RanecuEngine::flatArray(const int size, double* vect)
238 {
239  const int index = seq;
240  long seed1 = table[index][0];
241  long seed2 = table[index][1];
242  int k1, k2;
243  register int i;
244 
245  for (i=0; i<size; ++i)
246  {
247  k1 = (int)(seed1/ecuyer_b);
248  k2 = (int)(seed2/ecuyer_e);
249 
250  seed1 = ecuyer_a*(seed1-k1*ecuyer_b)-k1*ecuyer_c;
251  if (seed1 < 0) seed1 += shift1;
252  seed2 = ecuyer_d*(seed2-k2*ecuyer_e)-k2*ecuyer_f;
253  if (seed2 < 0) seed2 += shift2;
254 
255  long diff = seed1-seed2;
256  if (diff <= 0) diff += (shift1-1);
257 
258  vect[i] = (double)(diff*prec);
259  }
260  table[index][0] = seed1;
261  table[index][1] = seed2;
262 }
263 
264 RanecuEngine::operator unsigned int() {
265  const int index = seq;
266  long seed1 = table[index][0];
267  long seed2 = table[index][1];
268 
269  int k1 = (int)(seed1/ecuyer_b);
270  int k2 = (int)(seed2/ecuyer_e);
271 
272  seed1 = ecuyer_a*(seed1-k1*ecuyer_b)-k1*ecuyer_c;
273  if (seed1 < 0) seed1 += shift1;
274  seed2 = ecuyer_d*(seed2-k2*ecuyer_e)-k2*ecuyer_f;
275  if (seed2 < 0) seed2 += shift2;
276 
277  table[index][0] = seed1;
278  table[index][1] = seed2;
279  long diff = seed1-seed2;
280  if( diff <= 0 ) diff += (shift1-1);
281 
282  return ((diff << 1) | (seed1&1))& 0xffffffff;
283 }
284 
285 std::ostream & RanecuEngine::put( std::ostream& os ) const
286 {
287  char beginMarker[] = "RanecuEngine-begin";
288  os << beginMarker << "\nUvec\n";
289  std::vector<unsigned long> v = put();
290  for (unsigned int i=0; i<v.size(); ++i) {
291  os << v[i] << "\n";
292  }
293  return os;
294 #ifdef REMOVED
295  char endMarker[] = "RanecuEngine-end";
296  os << " " << beginMarker << "\n";
297  os << theSeed << " ";
298  for (int i=0; i<2; ++i) {
299  os << table[theSeed][i] << "\n";
300  }
301  os << endMarker << "\n";
302  return os;
303 #endif
304 }
305 
306 std::vector<unsigned long> RanecuEngine::put () const {
307  std::vector<unsigned long> v;
308  v.push_back (engineIDulong<RanecuEngine>());
309  v.push_back(static_cast<unsigned long>(theSeed));
310  v.push_back(static_cast<unsigned long>(table[theSeed][0]));
311  v.push_back(static_cast<unsigned long>(table[theSeed][1]));
312  return v;
313 }
314 
315 std::istream & RanecuEngine::get ( std::istream& is )
316 {
317  char beginMarker [MarkerLen];
318 
319  is >> std::ws;
320  is.width(MarkerLen); // causes the next read to the char* to be <=
321  // that many bytes, INCLUDING A TERMINATION \0
322  // (Stroustrup, section 21.3.2)
323  is >> beginMarker;
324  if (strcmp(beginMarker,"RanecuEngine-begin")) {
325  is.clear(std::ios::badbit | is.rdstate());
326  std::cerr << "\nInput stream mispositioned or"
327  << "\nRanecuEngine state description missing or"
328  << "\nwrong engine type found." << std::endl;
329  return is;
330  }
331  return getState(is);
332 }
333 
334 std::string RanecuEngine::beginTag ( ) {
335  return "RanecuEngine-begin";
336 }
337 
338 std::istream & RanecuEngine::getState ( std::istream& is )
339 {
340  if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) {
341  std::vector<unsigned long> v;
342  unsigned long uu;
343  for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
344  is >> uu;
345  if (!is) {
346  is.clear(std::ios::badbit | is.rdstate());
347  std::cerr << "\nRanecuEngine state (vector) description improper."
348  << "\ngetState() has failed."
349  << "\nInput stream is probably mispositioned now." << std::endl;
350  return is;
351  }
352  v.push_back(uu);
353  }
354  getState(v);
355  return (is);
356  }
357 
358 // is >> theSeed; Removed, encompassed by possibleKeywordInput()
359  char endMarker [MarkerLen];
360  for (int i=0; i<2; ++i) {
361  is >> table[theSeed][i];
362  }
363  is >> std::ws;
364  is.width(MarkerLen);
365  is >> endMarker;
366  if (strcmp(endMarker,"RanecuEngine-end")) {
367  is.clear(std::ios::badbit | is.rdstate());
368  std::cerr << "\nRanecuEngine state description incomplete."
369  << "\nInput stream is probably mispositioned now." << std::endl;
370  return is;
371  }
372 
373  seq = int(theSeed);
374  return is;
375 }
376 
377 bool RanecuEngine::get (const std::vector<unsigned long> & v) {
378  if ((v[0] & 0xffffffffUL) != engineIDulong<RanecuEngine>()) {
379  std::cerr <<
380  "\nRanecuEngine get:state vector has wrong ID word - state unchanged\n";
381  return false;
382  }
383  return getState(v);
384 }
385 
386 bool RanecuEngine::getState (const std::vector<unsigned long> & v) {
387  if (v.size() != VECTOR_STATE_SIZE ) {
388  std::cerr <<
389  "\nRanecuEngine get:state vector has wrong length - state unchanged\n";
390  return false;
391  }
392  theSeed = v[1];
393  table[theSeed][0] = v[2];
394  table[theSeed][1] = v[3];
395  seq = int(theSeed);
396  return true;
397 }
398 
399 
400 } // namespace CLHEP
double
#define double(obj)
Definition: excDblThrow.cc:32
CLHEP::RanecuEngine::setSeed
void setSeed(long index, int dum=0)
Definition: RanecuEngine.cc:106
CLHEP::HepRandomEngine
Definition: Matrix/CLHEP/Random/RandomEngine.h:55
CLHEP::RanecuEngine::flat
double flat()
Definition: RanecuEngine.cc:214
CLHEP::HepRandomEngine::theSeed
long theSeed
Definition: Matrix/CLHEP/Random/RandomEngine.h:144
CLHEP::RanecuEngine::ecuyer_a
static const int ecuyer_a
Definition: Matrix/CLHEP/Random/RanecuEngine.h:107
CLHEP::HepRandom::getTheTableSeeds
static void getTheTableSeeds(long *seeds, int index)
Definition: Random.cc:152
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::RanecuEngine::shift2
static const int shift2
Definition: Matrix/CLHEP/Random/RanecuEngine.h:114
CLHEP::RanecuEngine::~RanecuEngine
virtual ~RanecuEngine()
Definition: RanecuEngine.cc:104
CLHEP::HepRandomEngine::theSeeds
const long * theSeeds
Definition: Matrix/CLHEP/Random/RandomEngine.h:145
CLHEP::RanecuEngine::ecuyer_d
static const int ecuyer_d
Definition: Matrix/CLHEP/Random/RanecuEngine.h:110
CLHEP::RanecuEngine::engineName
static std::string engineName()
Definition: Matrix/CLHEP/Random/RanecuEngine.h:97
CLHEP::RanecuEngine::flatArray
void flatArray(const int size, double *vect)
Definition: RanecuEngine.cc:237
CLHEP::RanecuEngine::setIndex
void setIndex(long index)
Definition: RanecuEngine.cc:128
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::RanecuEngine::RanecuEngine
RanecuEngine()
Definition: RanecuEngine.cc:65
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::RanecuEngine::ecuyer_c
static const int ecuyer_c
Definition: Matrix/CLHEP/Random/RanecuEngine.h:109
CLHEP::RanecuEngine::saveStatus
void saveStatus(const char filename[]="Ranecu.conf") const
Definition: RanecuEngine.cc:135
CLHEP::RanecuEngine::restoreStatus
void restoreStatus(const char filename[]="Ranecu.conf")
Definition: RanecuEngine.cc:166
CLHEP::RanecuEngine::showStatus
void showStatus() const
Definition: RanecuEngine.cc:203
CLHEP::RanecuEngine::ecuyer_f
static const int ecuyer_f
Definition: Matrix/CLHEP/Random/RanecuEngine.h:112
CLHEP::possibleKeywordInput
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
Definition: Matrix/CLHEP/Random/RandomEngine.h:168
CLHEP::RanecuEngine::ecuyer_e
static const int ecuyer_e
Definition: Matrix/CLHEP/Random/RanecuEngine.h:111
j
long j
Definition: JamesRandomSeeding.txt:28
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
CLHEP::RanecuEngine::shift1
static const int shift1
Definition: Matrix/CLHEP/Random/RanecuEngine.h:113
CLHEP::RanecuEngine::getState
virtual std::istream & getState(std::istream &is)
Definition: RanecuEngine.cc:338
CLHEP::RanecuEngine::setSeeds
void setSeeds(const long *seeds, int index=-1)
Definition: RanecuEngine.cc:116
CLHEP::RanecuEngine::beginTag
static std::string beginTag()
Definition: RanecuEngine.cc:334
i
long i
Definition: JamesRandomSeeding.txt:27
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
CLHEP::RanecuEngine::name
std::string name() const
Definition: RanecuEngine.cc:54
CLHEP::RanecuEngine::VECTOR_STATE_SIZE
static const unsigned int VECTOR_STATE_SIZE
Definition: Matrix/CLHEP/Random/RanecuEngine.h:116
CLHEP::RanecuEngine::ecuyer_b
static const int ecuyer_b
Definition: Matrix/CLHEP/Random/RanecuEngine.h:108
CLHEP::RanecuEngine::put
std::vector< unsigned long > put() const
Definition: RanecuEngine.cc:306
CLHEP::RanecuEngine::get
virtual std::istream & get(std::istream &is)
Definition: RanecuEngine.cc:315