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