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

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