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

DRand48Engine.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 // $Id: DRand48Engine.cc,v 1.7 2010/07/29 16:50:34 garren Exp $
3 // -----------------------------------------------------------------------
4 // HEP Random
5 // --- DRand48Engine ---
6 // class implementation file
7 // -----------------------------------------------------------------------
8 // This file is part of Geant4 (simulation toolkit for HEP).
9 
10 // =======================================================================
11 // G.Cosmo - Created: 5th September 1995
12 // - Minor corrections: 31st October 1996
13 // - Added methods for engine status: 19th November 1996
14 // - Added srand48(), seed48(), drand48() implementations
15 // for Windows/NT: 6th March 1997
16 // - Fixed bug in setSeeds(): 15th September 1997
17 // - Private copy constructor and operator=: 26th Feb 1998
18 // J.Marraffino - Added stream operators and related constructor.
19 // Added automatic seed selection from seed table and
20 // engine counter: 16th Feb 1998
21 // J.Marraffino - Remove dependence on hepString class 13 May 1999
22 // E.Tcherniaev - More accurate code for drand48() on NT base on
23 // a code extracted from GNU C Library 2.1.3: 8th Nov 2000
24 // M. Fischler - In restore, checkFile for file not found 03 Dec 2004
25 // M. Fischler - Methods for distrib. instacne save/restore 12/8/04
26 // M. Fischler - split get() into tag validation and
27 // getState() for anonymous restores 12/27/04
28 // M. Fischler - put/get for vectors of ulongs 3/8/05
29 // M. Fischler - State-saving using only ints, for portability 4/12/05
30 //
31 // =======================================================================
32 
33 #include "CLHEP/Random/defs.h"
34 #include "CLHEP/Random/Random.h"
35 #include "CLHEP/Random/DRand48Engine.h"
36 #include "CLHEP/Random/RandomFunc.h"
37 #include "CLHEP/Random/engineIDulong.h"
38 #include <string.h> // for strcmp
39 #include <stdlib.h> // for std::abs(int)
40 
41 //#define TRACE_IO
42 
43 namespace CLHEP {
44 
45 static const int MarkerLen = 64; // Enough room to hold a begin or end marker.
46 // Number of instances with automatic seed selection
47 int DRand48Engine::numEngines = 0;
48 
49 std::string DRand48Engine::name() const {return "DRand48Engine";}
50 
51 // Maximum index into the seed table
52 int DRand48Engine::maxIndex = 215;
53 
56 {
57  setSeed(seed,0);
58  setSeeds(&theSeed,0);
59 }
60 
63 {
64  long seeds[2];
65  long seed;
66 
67  int cycle = abs(int(numEngines/maxIndex));
68  int curIndex = abs(int(numEngines%maxIndex));
69  ++numEngines;
70  long mask = ((cycle & 0x007fffff) << 8);
71  HepRandom::getTheTableSeeds( seeds, curIndex );
72  seed = seeds[0]^mask;
73  setSeed(seed,0);
74  setSeeds(&theSeed,0);
75 }
76 
77 DRand48Engine::DRand48Engine(int rowIndex, int colIndex)
79 {
80  long seed;
81  long seeds[2];
82 
83  int cycle = abs(int(rowIndex/maxIndex));
84  int row = abs(int(rowIndex%maxIndex));
85  int col = abs(int(colIndex%2));
86  long mask = ((cycle & 0x000007ff) << 20);
87  HepRandom::getTheTableSeeds( seeds, row );
88  seed = (seeds[col])^mask;
89  setSeed(seed,0);
90  setSeeds(&theSeed,0);
91 }
92 
95 {
96  is >> *this;
97 }
98 
100 
101 void DRand48Engine::setSeed(long seed, int)
102 {
103  srand48( seed );
104  theSeed = seed;
105 }
106 
107 void DRand48Engine::setSeeds(const long* seeds, int)
108 {
109  setSeed(seeds ? *seeds : 19780503L, 0);
110  theSeeds = seeds;
111 }
112 
113 void DRand48Engine::saveStatus( const char filename[] ) const
114 {
115  std::ofstream outFile( filename, std::ios::out ) ;
116 
117  if (!outFile.bad()) {
118  outFile << "Uvec\n";
119  std::vector<unsigned long> v = put();
120  #ifdef TRACE_IO
121  std::cout << "Result of v = put() is:\n";
122  #endif
123  for (unsigned int i=0; i<v.size(); ++i) {
124  outFile << v[i] << "\n";
125  #ifdef TRACE_IO
126  std::cout << v[i] << " ";
127  if (i%6==0) std::cout << "\n";
128  #endif
129  }
130  #ifdef TRACE_IO
131  std::cout << "\n";
132  #endif
133  }
134 
135 #ifdef REMOVED
136  unsigned short dummy[] = { 0, 0, 0 };
137  unsigned short* cseed = seed48(dummy);
138  if (!outFile.bad()) {
139  outFile << theSeed << std::endl;
140  for (int i=0; i<3; ++i) {
141  outFile << cseed[i] << std::endl;
142  dummy[i] = cseed[i];
143  }
144  seed48(dummy);
145  }
146 #endif
147 }
148 
149 void DRand48Engine::restoreStatus( const char filename[] )
150 {
151  std::ifstream inFile( filename, std::ios::in);
152  unsigned short cseed[3];
153 
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  #ifdef TRACE_IO
164  std::cout << "ivec = " << ivec << " xin = " << xin << " ";
165  if (ivec%3 == 0) std::cout << "\n";
166  #endif
167  if (!inFile) {
168  inFile.clear(std::ios::badbit | inFile.rdstate());
169  std::cerr << "\nDRand48Engine state (vector) description improper."
170  << "\nrestoreStatus has failed."
171  << "\nInput stream is probably mispositioned now." << std::endl;
172  return;
173  }
174  v.push_back(xin);
175  }
176  getState(v);
177  return;
178  }
179 
180  if (!inFile.bad() && !inFile.eof()) {
181  inFile >> theSeed;
182  for (int i=0; i<3; ++i)
183 // inFile >> theSeed; removed -- encompased by possibleKeywordInput
184  seed48(cseed);
185  }
186 }
187 
189 {
190  unsigned short dummy[] = { 0, 0, 0 };
191  unsigned short* cseed = seed48(dummy);
192  std::cout << std::endl;
193  std::cout << "-------- DRand48 engine status ---------" << std::endl;
194  std::cout << " Initial seed = " << theSeed << std::endl;
195  std::cout << " Current seeds = " << cseed[0] << ", ";
196  std::cout << cseed[1] << ", ";
197  std::cout << cseed[2] << std::endl;
198  std::cout << "----------------------------------------" << std::endl;
199  for (int i=0; i<3; ++i)
200  dummy[i] = cseed[i];
201  seed48(dummy);
202 }
203 
205 {
206  double num = 0.;
207 
208  while (num == 0.)
209  num = drand48();
210  return num;
211 }
212 
213 void DRand48Engine::flatArray(const int size, double* vect)
214 {
215  int i;
216 
217  for (i=0; i<size; ++i)
218  vect[i]=flat();
219 }
220 
221 std::ostream & DRand48Engine::put ( std::ostream& os ) const
222 {
223  char beginMarker[] = "DRand48Engine-begin";
224  os << beginMarker << "\nUvec\n";
225  std::vector<unsigned long> v = put();
226  for (unsigned int i=0; i<v.size(); ++i) {
227  os << v[i] << "\n";
228  }
229  return os;
230 
231 #ifdef REMOVED
232  unsigned short dummy[] = { 0, 0, 0 };
233  unsigned short* cseed = seed48(dummy);
234  char endMarker[] = "DRand48Engine-end";
235  os << " " << beginMarker << " ";
236  os << theSeed << " ";
237  for (int i=0; i<3; ++i) {
238  dummy[i] = cseed[i];
239  os << cseed[i] << " ";
240  }
241  os << endMarker << " ";
242  seed48(dummy);
243  return os;
244 #endif
245 }
246 
247 std::vector<unsigned long> DRand48Engine::put () const {
248  std::vector<unsigned long> v;
249  v.push_back (engineIDulong<DRand48Engine>());
250  unsigned short dummy[] = { 0, 0, 0 };
251  unsigned short* cseed = seed48(dummy);
252  for (int i=0; i<3; ++i) {
253  dummy[i] = cseed[i];
254  v.push_back (static_cast<unsigned long>(cseed[i]));
255  }
256  seed48(dummy);
257  return v;
258 }
259 
260 std::istream & DRand48Engine::get ( std::istream& is )
261 {
262  char beginMarker [MarkerLen];
263  is >> std::ws;
264  is.width(MarkerLen); // causes the next read to the char* to be <=
265  // that many bytes, INCLUDING A TERMINATION \0
266  // (Stroustrup, section 21.3.2)
267  is >> beginMarker;
268  if (strcmp(beginMarker,"DRand48Engine-begin")) {
269  is.clear(std::ios::badbit | is.rdstate());
270  std::cerr << "\nInput stream mispositioned or"
271  << "\nDRand48Engine state description missing or"
272  << "\nwrong engine type found." << std::endl;
273  return is;
274  }
275  return getState(is);
276 }
277 
278 std::string DRand48Engine::beginTag ( ) {
279  return "DRand48Engine-begin";
280 }
281 
282 std::istream & DRand48Engine::getState ( std::istream& is )
283 {
284  unsigned short cseed[3];
285  if ( possibleKeywordInput ( is, "Uvec", cseed[0] ) ) {
286  std::vector<unsigned long> v;
287  unsigned long uu;
288  #ifdef TRACE_IO
289  std::cout << "DRand48Engine::getState detected Uvec keyword\n";
290  uu = 999999;
291  #endif
292  for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
293  uu = 999999;
294  is >> uu;
295  #ifdef TRACE_IO
296  std::cout << "ivec = " << ivec << " uu = " << uu << "\n";
297  #endif
298  if (!is) {
299  is.clear(std::ios::badbit | is.rdstate());
300  std::cerr << "\nDRand48Engine state (vector) description improper."
301  << "\ngetState() has failed."
302  << "\nInput stream is probably mispositioned now." << std::endl;
303  return is;
304  }
305  v.push_back(uu);
306  }
307  getState(v);
308  return (is);
309  }
310 
311 // is >> cseed[0] was removed from loop, encompassed by possibleKeywordInput()
312 
313  char endMarker [MarkerLen];
314  is >> theSeed;
315  for (int i=1; i<3; ++i) {
316  is >> cseed[i];
317  }
318  is >> std::ws;
319  is.width(MarkerLen);
320  is >> endMarker;
321  if (strcmp(endMarker,"DRand48Engine-end")) {
322  is.clear(std::ios::badbit | is.rdstate());
323  std::cerr << "\nDRand48Engine state description incomplete."
324  << "\nInput stream is probably mispositioned now." << std::endl;
325  return is;
326  }
327  seed48(cseed);
328  return is;
329 }
330 
331 bool DRand48Engine::get (const std::vector<unsigned long> & v) {
332  if ((v[0] & 0xffffffffUL) != engineIDulong<DRand48Engine>()) {
333  std::cerr <<
334  "\nDRand48Engine get:state vector has wrong ID word - state unchanged\n";
335  return false;
336  }
337  return getState(v);
338 }
339 
340 bool DRand48Engine::getState (const std::vector<unsigned long> & v) {
341  if (v.size() != VECTOR_STATE_SIZE ) {
342  std::cerr <<
343  "\nDRand48Engine getState:state vector has wrong length - state unchanged\n";
344  return false;
345  }
346  unsigned short cseed[3];
347  for (int i=0; i<3; ++i) {
348  cseed[i] = static_cast<unsigned short>(v[i+1]);
349  }
350  seed48(cseed);
351  return true;
352 }
353 
354 } // namespace CLHEP