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