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);
79  HepRandom::getTheTableSeeds( seeds, curIndex );
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);
95  HepRandom::getTheTableSeeds( seeds, row );
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
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
void setSeed(long seed, int dum=0)
Definition: JamesRandom.cc:201
void saveStatus(const char filename[]="JamesRand.conf") const
Definition: JamesRandom.cc:109
static std::string beginTag()
Definition: JamesRandom.cc:362
void flatArray(const int size, double *vect)
Definition: JamesRandom.cc:284
std::vector< unsigned long > put() const
Definition: JamesRandom.cc:326
static double longs2double(const std::vector< unsigned long > &v)
Definition: DoubConv.cc:106
std::string name() const
Definition: JamesRandom.cc:54
void setSeeds(const long *seeds, int dum=0)
Definition: JamesRandom.cc:253
virtual std::istream & getState(std::istream &is)
Definition: JamesRandom.cc:366
static const unsigned int VECTOR_STATE_SIZE
static std::vector< unsigned long > dto2longs(double d)
Definition: DoubConv.cc:90
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
void showStatus() const
Definition: JamesRandom.cc:185
void restoreStatus(const char filename[]="JamesRand.conf")
Definition: JamesRandom.cc:143
virtual std::istream & get(std::istream &is)
Definition: JamesRandom.cc:345