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) {
75  HepRandom::getTheTableSeeds(table[j],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) {
90  HepRandom::getTheTableSeeds(table[j],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 
98 RanecuEngine::RanecuEngine(std::istream& is)
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
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
void restoreStatus(const char filename[]="Ranecu.conf")
#define double(obj)
Definition: excDblThrow.cc:32
void showStatus() const
void setSeeds(const long *seeds, int index=-1)
void setSeed(long index, int dum=0)
virtual std::istream & getState(std::istream &is)
void setIndex(long index)
static const unsigned int VECTOR_STATE_SIZE
virtual std::istream & get(std::istream &is)
std::string name() const
Definition: RanecuEngine.cc:54
void flatArray(const int size, double *vect)
static bool checkFile(std::istream &file, const std::string &filename, const std::string &classname, const std::string &methodname)
Definition: RandomEngine.cc:46
static void getTheTableSeeds(long *seeds, int index)
Definition: Random.cc:152
std::vector< unsigned long > put() const
void saveStatus(const char filename[]="Ranecu.conf") const
static std::string beginTag()