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

JamesRandom.cc
Go to the documentation of this file.
1 // $Id: JamesRandom.cc,v 1.6 2010/06/16 17:24:53 garren Exp $
2 // -*- C++ -*-
3 //
4 // -----------------------------------------------------------------------
5 // HEP Random
6 // --- HepJamesRandom ---
7 // class implementation file
8 // -----------------------------------------------------------------------
9 // This file is part of Geant4 (simulation toolkit for HEP).
10 //
11 // This algorithm implements the original universal random number generator
12 // as proposed by Marsaglia & Zaman in report FSU-SCRI-87-50 and coded
13 // in FORTRAN77 by Fred James as the RANMAR generator, part of the MATHLIB
14 // HEP library.
15 
16 // =======================================================================
17 // Gabriele Cosmo - Created: 5th September 1995
18 // - Fixed a bug in setSeed(): 26th February 1996
19 // - Minor corrections: 31st October 1996
20 // - Added methods for engine status: 19th November 1996
21 // - Fixed bug in setSeeds(): 15th September 1997
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 // V. Innocente - changed pointers to indices 3 may 2000
28 // M. Fischler - In restore, checkFile for file not found 03 Dec 2004
29 // M. Fischler - Methods for distrib. instacne 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 - Enforcement that seeds be non-negative
33 // (lest the sequence be non-random) 2/14/05
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/defs.h"
40 #include "CLHEP/Random/Random.h"
41 #include "CLHEP/Random/JamesRandom.h"
42 #include "CLHEP/Random/engineIDulong.h"
43 #include "CLHEP/Random/DoubConv.hh"
44 #include <string.h> // for strcmp
45 #include <cmath>
46 #include <cstdlib>
47 
48 //#define TRACE_IO
49 
50 namespace CLHEP {
51 
52 static const int MarkerLen = 64; // Enough room to hold a begin or end marker.
53 
54 std::string HepJamesRandom::name() const {return "HepJamesRandom";}
55 
56 // Number of instances with automatic seed selection
57 int HepJamesRandom::numEngines = 0;
58 
59 // Maximum index into the seed table
60 int HepJamesRandom::maxIndex = 215;
61 
64 {
65  setSeed(seed,0);
66  setSeeds(&theSeed,0);
67 }
68 
69 HepJamesRandom::HepJamesRandom() // 15 Feb. 1998 JMM
71 {
72  long seeds[2];
73  long seed;
74 
75  int cycle = std::abs(int(numEngines/maxIndex));
76  int curIndex = std::abs(int(numEngines%maxIndex));
77  ++numEngines;
78  long mask = ((cycle & 0x007fffff) << 8);
80  seed = seeds[0]^mask;
81  setSeed(seed,0);
82  setSeeds(&theSeed,0);
83 }
84 
85 HepJamesRandom::HepJamesRandom(int rowIndex, int colIndex) // 15 Feb. 1998 JMM
87 {
88  long seed;
89  long seeds[2];
90 
91  int cycle = std::abs(int(rowIndex/maxIndex));
92  int row = std::abs(int(rowIndex%maxIndex));
93  int col = std::abs(int(colIndex%2));
94  long mask = ((cycle & 0x000007ff) << 20);
96  seed = (seeds[col])^mask;
97  setSeed(seed,0);
98  setSeeds(&theSeed,0);
99 }
100 
102 : HepRandomEngine()
103 {
104  is >> *this;
105 }
106 
108 
109 void HepJamesRandom::saveStatus( const char filename[] ) const
110 {
111  std::ofstream outFile( filename, std::ios::out ) ;
112 
113  if (!outFile.bad()) {
114  outFile << "Uvec\n";
115  std::vector<unsigned long> v = put();
116  #ifdef TRACE_IO
117  std::cout << "Result of v = put() is:\n";
118  #endif
119  for (unsigned int i=0; i<v.size(); ++i) {
120  outFile << v[i] << "\n";
121  #ifdef TRACE_IO
122  std::cout << v[i] << " ";
123  if (i%6==0) std::cout << "\n";
124  #endif
125  }
126  #ifdef TRACE_IO
127  std::cout << "\n";
128  #endif
129  }
130 #ifdef REMOVED
131  int pos = j97;
132  outFile << theSeed << std::endl;
133  for (int i=0; i<97; ++i)
134  outFile << std::setprecision(20) << u[i] << " ";
135  outFile << std::endl;
136  outFile << std::setprecision(20) << c << " ";
137  outFile << std::setprecision(20) << cd << " ";
138  outFile << std::setprecision(20) << cm << std::endl;
139  outFile << pos << std::endl;
140 #endif
141 }
142 
143 void HepJamesRandom::restoreStatus( const char filename[] )
144 {
145  int ipos, jpos;
146  std::ifstream inFile( filename, std::ios::in);
147  if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) {
148  std::cerr << " -- Engine state remains unchanged\n";
149  return;
150  }
151  if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) {
152  std::vector<unsigned long> v;
153  unsigned long xin;
154  for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
155  inFile >> xin;
156  #ifdef TRACE_IO
157  std::cout << "ivec = " << ivec << " xin = " << xin << " ";
158  if (ivec%3 == 0) std::cout << "\n";
159  #endif
160  if (!inFile) {
161  inFile.clear(std::ios::badbit | inFile.rdstate());
162  std::cerr << "\nJamesRandom state (vector) description improper."
163  << "\nrestoreStatus has failed."
164  << "\nInput stream is probably mispositioned now." << std::endl;
165  return;
166  }
167  v.push_back(xin);
168  }
169  getState(v);
170  return;
171  }
172 
173  if (!inFile.bad() && !inFile.eof()) {
174 // inFile >> theSeed; removed -- encompased by possibleKeywordInput
175  for (int i=0; i<97; ++i)
176  inFile >> u[i];
177  inFile >> c; inFile >> cd; inFile >> cm;
178  inFile >> jpos;
179  ipos = (64+jpos)%97;
180  i97 = ipos;
181  j97 = jpos;
182  }
183 }
184 
186 {
187  std::cout << std::endl;
188  std::cout << "----- HepJamesRandom engine status -----" << std::endl;
189  std::cout << " Initial seed = " << theSeed << std::endl;
190  std::cout << " u[] = ";
191  for (int i=0; i<97; ++i)
192  std::cout << u[i] << " ";
193  std::cout << std::endl;
194  std::cout << " c = " << c << ", cd = " << cd << ", cm = " << cm
195  << std::endl;
196  std::cout << " i97 = " << i97 << ", u[i97] = " << u[i97] << std::endl;
197  std::cout << " j97 = " << j97 << ", u[j97] = " << u[j97] << std::endl;
198  std::cout << "----------------------------------------" << std::endl;
199 }
200 
201 void HepJamesRandom::setSeed(long seed, int)
202 {
203  // The input value for "seed" should be within the range [0,900000000]
204  //
205  // Negative seeds result in serious flaws in the randomness;
206  // seeds above 900000000 are OK because of the %177 in the expression for i,
207  // but may have the same effect as other seeds below 900000000.
208 
209  int m, n;
210  float s, t;
211  long mm;
212 
213  if (seed < 0) {
214  std::cout << "Seed for HepJamesRandom must be non-negative\n"
215  << "Seed value supplied was " << seed
216  << "\nUsing its absolute value instead\n";
217  seed = -seed;
218  }
219 
220  long ij = seed/30082;
221  long kl = seed - 30082*ij;
222  long i = (ij/177) % 177 + 2;
223  long j = ij % 177 + 2;
224  long k = (kl/169) % 178 + 1;
225  long l = kl % 169;
226 
227  theSeed = seed;
228 
229  for ( n = 1 ; n < 98 ; n++ ) {
230  s = 0.0;
231  t = 0.5;
232  for ( m = 1 ; m < 25 ; m++) {
233  mm = ( ( (i*j) % 179 ) * k ) % 179;
234  i = j;
235  j = k;
236  k = mm;
237  l = ( 53 * l + 1 ) % 169;
238  if ( (l*mm % 64 ) >= 32 )
239  s += t;
240  t *= 0.5;
241  }
242  u[n-1] = s;
243  }
244  c = 362436.0 / 16777216.0;
245  cd = 7654321.0 / 16777216.0;
246  cm = 16777213.0 / 16777216.0;
247 
248  i97 = 96;
249  j97 = 32;
250 
251 }
252 
253 void HepJamesRandom::setSeeds(const long* seeds, int)
254 {
255  setSeed(seeds ? *seeds : 19780503L, 0);
256  theSeeds = seeds;
257 }
258 
260 {
261  double uni;
262 
263  do {
264  uni = u[i97] - u[j97];
265  if ( uni < 0.0 ) uni++;
266  u[i97] = uni;
267 
268  if (i97 == 0) i97 = 96;
269  else i97--;
270 
271  if (j97 == 0) j97 = 96;
272  else j97--;
273 
274  c -= cd;
275  if (c < 0.0) c += cm;
276 
277  uni -= c;
278  if (uni < 0.0) uni += 1.0;
279  } while ( uni <= 0.0 || uni >= 1.0 );
280 
281  return uni;
282 }
283 
284 void HepJamesRandom::flatArray(const int size, double* vect)
285 {
286 // double uni;
287  int i;
288 
289  for (i=0; i<size; ++i) {
290  vect[i] = flat();
291  }
292 }
293 
294 HepJamesRandom::operator unsigned int() {
295  return ((unsigned int)(flat() * exponent_bit_32()) & 0xffffffff ) |
296  (((unsigned int)( u[i97] * exponent_bit_32())>>16) & 0xff);
297 }
298 
299 std::ostream & HepJamesRandom::put ( std::ostream& os ) const {
300  char beginMarker[] = "JamesRandom-begin";
301  os << beginMarker << "\nUvec\n";
302  std::vector<unsigned long> v = put();
303  for (unsigned int i=0; i<v.size(); ++i) {
304  os << v[i] << "\n";
305  }
306  return os;
307 #ifdef REMOVED
308  char endMarker[] = "JamesRandom-end";
309  int pos = j97;
310  int pr = os.precision(20);
311  os << " " << beginMarker << " ";
312  os << theSeed << " ";
313  for (int i=0; i<97; ++i) {
314  os << std::setprecision(20) << u[i] << "\n";
315  }
316  os << std::setprecision(20) << c << " ";
317  os << std::setprecision(20) << cd << " ";
318  os << std::setprecision(20) << cm << " ";
319  os << pos << "\n";
320  os << endMarker << "\n";
321  os.precision(pr);
322  return os;
323 #endif
324 }
325 
326 std::vector<unsigned long> HepJamesRandom::put () const {
327  std::vector<unsigned long> v;
328  v.push_back (engineIDulong<HepJamesRandom>());
329  std::vector<unsigned long> t;
330  for (int i=0; i<97; ++i) {
331  t = DoubConv::dto2longs(u[i]);
332  v.push_back(t[0]); v.push_back(t[1]);
333  }
334  t = DoubConv::dto2longs(c);
335  v.push_back(t[0]); v.push_back(t[1]);
336  t = DoubConv::dto2longs(cd);
337  v.push_back(t[0]); v.push_back(t[1]);
338  t = DoubConv::dto2longs(cm);
339  v.push_back(t[0]); v.push_back(t[1]);
340  v.push_back(static_cast<unsigned long>(j97));
341  return v;
342 }
343 
344 
345 std::istream & HepJamesRandom::get ( std::istream& is) {
346  char beginMarker [MarkerLen];
347  is >> std::ws;
348  is.width(MarkerLen); // causes the next read to the char* to be <=
349  // that many bytes, INCLUDING A TERMINATION \0
350  // (Stroustrup, section 21.3.2)
351  is >> beginMarker;
352  if (strcmp(beginMarker,"JamesRandom-begin")) {
353  is.clear(std::ios::badbit | is.rdstate());
354  std::cerr << "\nInput stream mispositioned or"
355  << "\nJamesRandom state description missing or"
356  << "\nwrong engine type found." << std::endl;
357  return is;
358  }
359  return getState(is);
360 }
361 
362 std::string HepJamesRandom::beginTag ( ) {
363  return "JamesRandom-begin";
364 }
365 
366 std::istream & HepJamesRandom::getState ( std::istream& is) {
367  if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) {
368  std::vector<unsigned long> v;
369  unsigned long uu;
370  for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
371  is >> uu;
372  if (!is) {
373  is.clear(std::ios::badbit | is.rdstate());
374  std::cerr << "\nJamesRandom state (vector) description improper."
375  << "\ngetState() has failed."
376  << "\nInput stream is probably mispositioned now." << std::endl;
377  return is;
378  }
379  v.push_back(uu);
380  }
381  getState(v);
382  return (is);
383  }
384 
385 // is >> theSeed; Removed, encompassed by possibleKeywordInput()
386 
387  int ipos, jpos;
388  char endMarker [MarkerLen];
389  for (int i=0; i<97; ++i) {
390  is >> u[i];
391  }
392  is >> c; is >> cd; is >> cm;
393  is >> jpos;
394  is >> std::ws;
395  is.width(MarkerLen);
396  is >> endMarker;
397  if(strcmp(endMarker,"JamesRandom-end")) {
398  is.clear(std::ios::badbit | is.rdstate());
399  std::cerr << "\nJamesRandom state description incomplete."
400  << "\nInput stream is probably mispositioned now." << std::endl;
401  return is;
402  }
403 
404  ipos = (64+jpos)%97;
405  i97 = ipos;
406  j97 = jpos;
407  return is;
408 }
409 
410 bool HepJamesRandom::get (const std::vector<unsigned long> & v) {
411  if ( (v[0] & 0xffffffffUL) != engineIDulong<HepJamesRandom>()) {
412  std::cerr <<
413  "\nHepJamesRandom get:state vector has wrong ID word - state unchanged\n";
414  return false;
415  }
416  return getState(v);
417 }
418 
419 bool HepJamesRandom::getState (const std::vector<unsigned long> & v) {
420  if (v.size() != VECTOR_STATE_SIZE ) {
421  std::cerr <<
422  "\nHepJamesRandom get:state vector has wrong length - state unchanged\n";
423  return false;
424  }
425  std::vector<unsigned long> t(2);
426  for (int i=0; i<97; ++i) {
427  t[0] = v[2*i+1]; t[1] = v[2*i+2];
428  u[i] = DoubConv::longs2double(t);
429  }
430  t[0] = v[195]; t[1] = v[196]; c = DoubConv::longs2double(t);
431  t[0] = v[197]; t[1] = v[198]; cd = DoubConv::longs2double(t);
432  t[0] = v[199]; t[1] = v[200]; cm = DoubConv::longs2double(t);
433  j97 = v[201];
434  i97 = (64+j97)%97;
435  return true;
436 }
437 
438 } // namespace CLHEP
CLHEP::HepJamesRandom::flat
double flat()
Definition: JamesRandom.cc:259
l
long l
Definition: JamesRandomSeeding.txt:30
CLHEP::HepRandomEngine
Definition: Matrix/CLHEP/Random/RandomEngine.h:55
CLHEP::HepRandomEngine::theSeed
long theSeed
Definition: Matrix/CLHEP/Random/RandomEngine.h:144
CLHEP::HepJamesRandom::~HepJamesRandom
virtual ~HepJamesRandom()
Definition: JamesRandom.cc:107
CLHEP::HepRandom::getTheTableSeeds
static void getTheTableSeeds(long *seeds, int index)
Definition: Random.cc:152
CLHEP::HepJamesRandom::put
std::vector< unsigned long > put() const
Definition: JamesRandom.cc:326
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::HepJamesRandom::setSeeds
void setSeeds(const long *seeds, int dum=0)
Definition: JamesRandom.cc:253
CLHEP::DoubConv::longs2double
static double longs2double(const std::vector< unsigned long > &v)
Definition: DoubConv.cc:106
CLHEP::HepJamesRandom::flatArray
void flatArray(const int size, double *vect)
Definition: JamesRandom.cc:284
CLHEP::detail::n
n
Definition: Ranlux64Engine.cc:85
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::HepJamesRandom::getState
virtual std::istream & getState(std::istream &is)
Definition: JamesRandom.cc:366
CLHEP::HepJamesRandom::HepJamesRandom
HepJamesRandom()
Definition: JamesRandom.cc:69
CLHEP::HepJamesRandom::name
std::string name() const
Definition: JamesRandom.cc:54
CLHEP
Definition: ClhepVersion.h:13
ij
ij
Definition: JamesRandomSeeding.txt:50
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::possibleKeywordInput
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
Definition: Matrix/CLHEP/Random/RandomEngine.h:168
CLHEP::HepJamesRandom::restoreStatus
void restoreStatus(const char filename[]="JamesRand.conf")
Definition: JamesRandom.cc:143
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::HepJamesRandom::get
virtual std::istream & get(std::istream &is)
Definition: JamesRandom.cc:345
CLHEP::HepJamesRandom::saveStatus
void saveStatus(const char filename[]="JamesRand.conf") const
Definition: JamesRandom.cc:109
s
Methods applicble to containers of as in std::list< LorentzVector > s
Definition: keyMergeIssues.doc:328
i
long i
Definition: JamesRandomSeeding.txt:27
CLHEP::DoubConv::dto2longs
static std::vector< unsigned long > dto2longs(double d)
Definition: DoubConv.cc:90
CLHEP::HepJamesRandom::setSeed
void setSeed(long seed, int dum=0)
Definition: JamesRandom.cc:201
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::HepJamesRandom::beginTag
static std::string beginTag()
Definition: JamesRandom.cc:362
CLHEP::HepJamesRandom::VECTOR_STATE_SIZE
static const unsigned int VECTOR_STATE_SIZE
Definition: Matrix/CLHEP/Random/JamesRandom.h:95
kl
long kl
Definition: JamesRandomSeeding.txt:26
CLHEP::HepJamesRandom::showStatus
void showStatus() const
Definition: JamesRandom.cc:185
k
long k
Definition: JamesRandomSeeding.txt:29
CLHEP::HepJamesRandom::engineName
static std::string engineName()
Definition: Matrix/CLHEP/Random/JamesRandom.h:89