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

Hurd288Engine.cc
Go to the documentation of this file.
1 // $Id: Hurd288Engine.cc,v 1.7 2010/07/20 18:07:17 garren Exp $
2 // -*- C++ -*-
3 //
4 // -----------------------------------------------------------------------
5 // HEP Random
6 // --- Hurd288Engine ---
7 // class implementation file
8 // -----------------------------------------------------------------------
9 // An interconnected shift register based on the paper presented by Hurd in
10 // IEEE Transactions on Computers c23, 2 Feb 1974.
11 // =======================================================================
12 // Ken Smith - Initial draft started: 23rd Jul 1998
13 // - Added conversion operators: 6th Aug 1998
14 // J. Marraffino - Added some explicit casts to deal with
15 // machines where sizeof(int) != sizeof(long) 22 Aug 1998
16 // M. Fischler - Modified use of the various exponents of 2
17 // to avoid per-instance space overhead and
18 // correct the rounding procedure 15 Sep 1998
19 // - Modified various methods to avoid any
20 // possibility of predicting the next number
21 // based on the last several: We skip one
22 // 32-bit word per cycle. 15 Sep 1998
23 // - modify word[0] in two of the constructors
24 // so that no sequence can ever accidentally
25 // be produced by differnt seeds. 15 Sep 1998
26 // J. Marraffino - Remove dependence on hepString class 13 May 1999
27 // M. Fischler - Put endl at end of a save 10 Apr 2001
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 - put/get for vectors of ulongs 3/14/05
33 // M. Fischler - State-saving using only ints, for portability 4/12/05
34 //
35 // =======================================================================
36 
37 #include "CLHEP/Random/defs.h"
38 #include "CLHEP/Random/Random.h"
39 #include "CLHEP/Random/Hurd288Engine.h"
40 #include "CLHEP/Random/engineIDulong.h"
41 #include <string.h> // for strcmp
42 #include <cstdlib> // for std::abs(int)
43 
44 using namespace std;
45 
46 namespace CLHEP {
47 
48 static const int MarkerLen = 64; // Enough room to hold a begin or end marker.
49 
50 std::string Hurd288Engine::name() const {return "Hurd288Engine";}
51 
52 // Number of instances with automatic seed selection
53 int Hurd288Engine::numEngines = 0;
54 
55 // Maximum index into the seed table
56 int Hurd288Engine::maxIndex = 215;
57 
58 static inline unsigned int f288(unsigned int a, unsigned int b, unsigned int c)
59 {
60  return ( ((b<<2) & 0x7ffc) | ((a<<2) & ~0x7ffc) | (a>>30) ) ^
61  ( (c<<1)|(c>>31) );
62 }
63 
64 Hurd288Engine::Hurd288Engine()
66 {
67  int cycle = std::abs(int(numEngines/maxIndex));
68  int curIndex = std::abs(int(numEngines%maxIndex));
69  long mask = ((cycle & 0x007fffff) << 8);
70  long seedlist[2];
71  HepRandom::getTheTableSeeds( seedlist, curIndex );
72  seedlist[0] ^= mask;
73  seedlist[1] = 0;
74  setSeeds(seedlist, 0);
75  words[0] ^= 0x1324abcd; // To make unique vs long or two unsigned
76  if (words[0]==0) words[0] = 1; // ints in the constructor
77  ++numEngines;
78  for( int i=0; i < 100; ++i ) flat(); // warm up just a bit
79 }
80 
81 Hurd288Engine::Hurd288Engine( std::istream& is )
83 {
84  is >> *this;
85 }
86 
89 {
90  long seedlist[2]={seed,0};
91  setSeeds(seedlist, 0);
92  words[0] ^= 0xa5482134; // To make unique vs default two unsigned
93  if (words[0]==0) words[0] = 1; // ints in the constructor
94  for( int i=0; i < 100; ++i ) flat(); // warm up just a bit
95 }
96 
97 Hurd288Engine::Hurd288Engine( int rowIndex, int colIndex )
99 {
100  int cycle = std::abs(int(rowIndex/maxIndex));
101  int row = std::abs(int(rowIndex%maxIndex));
102  int col = colIndex & 0x1;
103  long mask = (( cycle & 0x000007ff ) << 20 );
104  long seedlist[2];
105  HepRandom::getTheTableSeeds( seedlist, row );
106  seedlist[0] = (seedlist[col])^mask;
107  seedlist[1]= 0;
108  setSeeds(seedlist, 0);
109  for( int i=0; i < 100; ++i ) flat(); // warm up just a bit
110 }
111 
113 
114 void Hurd288Engine::advance() {
115 
116  register unsigned int W0, W1, W2, W3, W4, W5, W6, W7, W8;
117 
118  W0 = words[0];
119  W1 = words[1];
120  W2 = words[2];
121  W3 = words[3];
122  W4 = words[4];
123  W5 = words[5];
124  W6 = words[6];
125  W7 = words[7];
126  W8 = words[8];
127  W1 ^= W0; W0 = f288(W2, W3, W0);
128  W2 ^= W1; W1 = f288(W3, W4, W1);
129  W3 ^= W2; W2 = f288(W4, W5, W2);
130  W4 ^= W3; W3 = f288(W5, W6, W3);
131  W5 ^= W4; W4 = f288(W6, W7, W4);
132  W6 ^= W5; W5 = f288(W7, W8, W5);
133  W7 ^= W6; W6 = f288(W8, W0, W6);
134  W8 ^= W7; W7 = f288(W0, W1, W7);
135  W0 ^= W8; W8 = f288(W1, W2, W8);
136  words[0] = W0 & 0xffffffff;
137  words[1] = W1 & 0xffffffff;
138  words[2] = W2 & 0xffffffff;
139  words[3] = W3 & 0xffffffff;
140  words[4] = W4 & 0xffffffff;
141  words[5] = W5 & 0xffffffff;
142  words[6] = W6 & 0xffffffff;
143  words[7] = W7 & 0xffffffff;
144  words[8] = W8 & 0xffffffff;
145  wordIndex = 9;
146 
147 } // advance()
148 
149 
151 
152  if( wordIndex <= 2 ) { // MF 9/15/98:
153  // skip word 0 and use two words per flat
154  advance();
155  }
156 
157  // LG 6/30/2010
158  // define the order of execution for --wordIndex
159  double x = words[--wordIndex] * twoToMinus_32() ; // most significant part
160  double y = (words[--wordIndex]>>11) * twoToMinus_53() + // fill in rest of bits
161  nearlyTwoToMinus_54(); // make sure non-zero
162  return x + y;
163 }
164 
165 void Hurd288Engine::flatArray( const int size, double* vect ) {
166  for (int i = 0; i < size; ++i) {
167  vect[i] = flat();
168  }
169 }
170 
171 void Hurd288Engine::setSeed( long seed, int ) {
172  words[0] = (unsigned int)seed;
173  for (wordIndex = 1; wordIndex < 9; ++wordIndex) {
174  words[wordIndex] = 69607 * words[wordIndex-1] + 54329;
175  }
176 }
177 
178 void Hurd288Engine::setSeeds( const long* seeds, int ) {
179  setSeed( *seeds ? *seeds : 32767, 0 );
180  theSeeds = seeds;
181 }
182 
183 void Hurd288Engine::saveStatus( const char filename[] ) const {
184  std::ofstream outFile(filename, std::ios::out);
185  if( !outFile.bad() ) {
186  outFile << "Uvec\n";
187  std::vector<unsigned long> v = put();
188  #ifdef TRACE_IO
189  std::cout << "Result of v = put() is:\n";
190  #endif
191  for (unsigned int i=0; i<v.size(); ++i) {
192  outFile << v[i] << "\n";
193  #ifdef TRACE_IO
194  std::cout << v[i] << " ";
195  if (i%6==0) std::cout << "\n";
196  #endif
197  }
198  #ifdef TRACE_IO
199  std::cout << "\n";
200  #endif
201  }
202 #ifdef REMOVED
203  outFile << std::setprecision(20) << theSeed << " ";
204  outFile << wordIndex << " ";
205  for( int i = 0; i < 9; ++i ) {
206  outFile << words[i] << " ";
207  }
208  outFile << std::endl;
209 #endif
210 }
211 
212 void Hurd288Engine::restoreStatus( const char filename[] ) {
213  std::ifstream inFile(filename, std::ios::in);
214  if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) {
215  std::cerr << " -- Engine state remains unchanged\n";
216  return;
217  }
218  if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) {
219  std::vector<unsigned long> v;
220  unsigned long xin;
221  for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
222  inFile >> xin;
223  #ifdef TRACE_IO
224  std::cout << "ivec = " << ivec << " xin = " << xin << " ";
225  if (ivec%3 == 0) std::cout << "\n";
226  #endif
227  if (!inFile) {
228  inFile.clear(std::ios::badbit | inFile.rdstate());
229  std::cerr << "\nHurd288Engine state (vector) description improper."
230  << "\nrestoreStatus has failed."
231  << "\nInput stream is probably mispositioned now." << std::endl;
232  return;
233  }
234  v.push_back(xin);
235  }
236  getState(v);
237  return;
238  }
239 
240  if( !inFile.bad() ) {
241 // inFile >> theSeed; removed -- encompased by possibleKeywordInput
242  inFile >> wordIndex;
243  for( int i = 0; i < 9; ++i ) {
244  inFile >> words[i];
245  }
246  }
247 }
248 
250  std::cout << std::setprecision(20) << std::endl;
251  std::cout << "----------- Hurd2 engine status ----------" << std::endl;
252  std::cout << "Initial seed = " << theSeed << std::endl;
253  std::cout << "Current index = " << wordIndex << std::endl;
254  std::cout << "Current words = " << std::endl;
255  for( int i = 0; i < 9 ; ++i ) {
256  std::cout << " " << words[i] << std::endl;
257  }
258  std::cout << "-------------------------------------------" << std::endl;
259 }
260 
261 Hurd288Engine::operator float() {
262  if( wordIndex <= 1 ) { // MF 9/15/98: skip word 0
263  advance();
264  }
265  return words[--wordIndex ] * twoToMinus_32();
266 }
267 
268 Hurd288Engine::operator unsigned int() {
269  if( wordIndex <= 1 ) { // MF 9/15/98: skip word 0
270  advance();
271  }
272  return words[--wordIndex];
273 }
274 
275 std::ostream& Hurd288Engine::put(std::ostream& os) const {
276  char beginMarker[] = "Hurd288Engine-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 #ifdef REMOVED
284  char endMarker[] = "Hurd288Engine-end";
285  int pr = os.precision(20);
286  os << " " << beginMarker << " ";
287  os << theSeed << " ";
288  os << wordIndex << " ";
289  for (int i = 0; i < 9; ++i) {
290  os << words[i] << "\n";
291  }
292  os << endMarker << "\n ";
293  os.precision(pr);
294  return os;
295 #endif
296 }
297 
298 std::vector<unsigned long> Hurd288Engine::put () const {
299  std::vector<unsigned long> v;
300  v.push_back (engineIDulong<Hurd288Engine>());
301  v.push_back(static_cast<unsigned long>(wordIndex));
302  for (int i = 0; i < 9; ++i) {
303  v.push_back(static_cast<unsigned long>(words[i]));
304  }
305  return v;
306 }
307 
308 
309 std::istream& Hurd288Engine::get(std::istream& is) {
310  char beginMarker [MarkerLen];
311  is >> std::ws;
312  is.width(MarkerLen); // causes the next read to the char* to be <=
313  // that many bytes, INCLUDING A TERMINATION \0
314  // (Stroustrup, section 21.3.2)
315  is >> beginMarker;
316  if (strcmp(beginMarker,"Hurd288Engine-begin")) {
317  is.clear(std::ios::badbit | is.rdstate());
318  std::cerr << "\nInput mispositioned or"
319  << "\nHurd288Engine state description missing or"
320  << "\nwrong engine type found." << std::endl;
321  return is;
322  }
323  return getState(is);
324 }
325 
326 std::string Hurd288Engine::beginTag ( ) {
327  return "Hurd288Engine-begin";
328 }
329 
330 std::istream& Hurd288Engine::getState(std::istream& is) {
331  if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) {
332  std::vector<unsigned long> v;
333  unsigned long uu;
334  for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
335  is >> uu;
336  if (!is) {
337  is.clear(std::ios::badbit | is.rdstate());
338  std::cerr << "\nHurd288Engine state (vector) description improper."
339  << "\ngetState() has failed."
340  << "\nInput stream is probably mispositioned now." << std::endl;
341  return is;
342  }
343  v.push_back(uu);
344  }
345  getState(v);
346  return (is);
347  }
348 
349 // is >> theSeed; Removed, encompassed by possibleKeywordInput()
350 
351  char endMarker [MarkerLen];
352  is >> wordIndex;
353  for (int i = 0; i < 9; ++i) {
354  is >> words[i];
355  }
356  is >> std::ws;
357  is.width(MarkerLen);
358  is >> endMarker;
359  if (strcmp(endMarker,"Hurd288Engine-end")) {
360  is.clear(std::ios::badbit | is.rdstate());
361  std::cerr << "\nHurd288Engine state description incomplete."
362  << "\nInput stream is probably mispositioned now." << std::endl;
363  return is;
364  }
365  return is;
366 }
367 
368 bool Hurd288Engine::get (const std::vector<unsigned long> & v) {
369  if ((v[0] & 0xffffffffUL) != engineIDulong<Hurd288Engine>()) {
370  std::cerr <<
371  "\nHurd288Engine get:state vector has wrong ID word - state unchanged\n";
372  std::cerr << "The correct ID would be " << engineIDulong<Hurd288Engine>()
373  << "; the actual ID is " << v[0] << "\n";
374  return false;
375  }
376  return getState(v);
377 }
378 
379 bool Hurd288Engine::getState (const std::vector<unsigned long> & v) {
380  if (v.size() != VECTOR_STATE_SIZE ) {
381  std::cerr <<
382  "\nHurd288Engine get:state vector has wrong length - state unchanged\n";
383  return false;
384  }
385  wordIndex = v[1];
386  for (int i = 0; i < 9; ++i) {
387  words[i] = v[i+2];
388  }
389  return true;
390 }
391 
392 } // namespace CLHEP