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

DualRand.cc
Go to the documentation of this file.
1 // $Id: DualRand.cc,v 1.5 2010/06/16 17:24:53 garren Exp $
2 // -*- C++ -*-
3 //
4 // -----------------------------------------------------------------------
5 // Hep Random
6 // --- DualRand ---
7 // class implementation file
8 // -----------------------------------------------------------------------
9 // Exclusive or of a feedback shift register and integer congruence
10 // random number generator. The feedback shift register uses offsets
11 // 127 and 97. The integer congruence generator uses a different
12 // multiplier for each stream. The multipliers are chosen to give
13 // full period and maximum "potency" for modulo 2^32. The period of
14 // the combined random number generator is 2^159 - 2^32, and the
15 // sequences are different for each stream (not just started in a
16 // different place).
17 //
18 // In the original generator used on ACPMAPS:
19 // The feedback shift register generates 24 bits at a time, and
20 // the high order 24 bits of the integer congruence generator are
21 // used.
22 //
23 // Instead, to conform with more modern engine concepts, we generate
24 // 32 bits at a time and use the full 32 bits of the congruence
25 // generator.
26 //
27 // References:
28 // Knuth
29 // Tausworthe
30 // Golomb
31 //=========================================================================
32 // Ken Smith - Removed std::pow() from flat() method: 21 Jul 1998
33 // - Added conversion operators: 6 Aug 1998
34 // J. Marraffino - Added some explicit casts to deal with
35 // machines where sizeof(int) != sizeof(long) 22 Aug 1998
36 // M. Fischler - Modified constructors taking seeds to not
37 // depend on numEngines (same seeds should
38 // produce same sequences). Default still
39 // depends on numEngines. 14 Sep 1998
40 // - Modified use of the various exponents of 2
41 // to avoid per-instance space overhead and
42 // correct the rounding procedure 15 Sep 1998
43 // J. Marraffino - Remove dependence on hepString class 13 May 1999
44 // M. Fischler - Put endl at end of a save 10 Apr 2001
45 // M. Fischler - In restore, checkFile for file not found 03 Dec 2004
46 // M. Fischler - methods for distrib. instacne save/restore 12/8/04
47 // M. Fischler - split get() into tag validation and
48 // getState() for anonymous restores 12/27/04
49 // Mark Fischler - methods for vector save/restore 3/7/05
50 // M. Fischler - State-saving using only ints, for portability 4/12/05
51 //
52 //=========================================================================
53 
54 #include "CLHEP/Random/DualRand.h"
55 #include "CLHEP/Random/defs.h"
56 #include "CLHEP/Random/engineIDulong.h"
57 #include <string.h> // for strcmp
58 
59 namespace CLHEP {
60 
61 static const int MarkerLen = 64; // Enough room to hold a begin or end marker.
62 
63 std::string DualRand::name() const {return "DualRand";}
64 
65 // Number of instances with automatic seed selection
66 int DualRand::numEngines = 0;
67 
68 // The following constructors (excluding the istream constructor) fill
69 // the bits of the tausworthe and the starting state of the integer
70 // congruence based on the seed. In addition, it sets up the multiplier
71 // for the integer congruence based on the stream number, so you have
72 // absolutely independent streams.
73 
75 : HepRandomEngine(),
76  tausworthe (1234567 + numEngines + 175321),
77  integerCong(69607 * tausworthe + 54329, numEngines)
78 {
79  theSeed = 1234567;
80  ++numEngines;
81 }
82 
84 : HepRandomEngine(),
85  tausworthe ((unsigned int)seed + 175321),
86  integerCong(69607 * tausworthe + 54329, 8043) // MF - not numEngines
87 {
88  theSeed = seed;
89 }
90 
91 DualRand::DualRand(std::istream & is)
93 {
94  is >> *this;
95 }
96 
97 DualRand::DualRand(int rowIndex, int colIndex)
98 : HepRandomEngine(),
99  tausworthe (rowIndex + 1000 * colIndex + 85329),
100  integerCong(69607 * tausworthe + 54329, 1123) // MF - not numengines
101 {
102  theSeed = rowIndex;
103 }
104 
106 
107 double DualRand::flat() {
108  unsigned int ic ( integerCong );
109  unsigned int t ( tausworthe );
110  return ( (t ^ ic) * twoToMinus_32() + // most significant part
111  (t >> 11) * twoToMinus_53() + // fill in remaining bits
112  nearlyTwoToMinus_54() // make sure non-zero
113  );
114 }
115 
116 void DualRand::flatArray(const int size, double* vect) {
117  for (int i = 0; i < size; ++i) {
118  vect[i] = flat();
119  }
120 }
121 
122 void DualRand::setSeed(long seed, int) {
123  theSeed = seed;
124  tausworthe = Tausworthe((unsigned int)seed + numEngines + 175321);
125  integerCong = IntegerCong(69607 * tausworthe + 54329, numEngines);
126 }
127 
128 void DualRand::setSeeds(const long * seeds, int) {
129  setSeed(seeds ? *seeds : 1234567, 0);
130  theSeeds = seeds;
131 }
132 
133 void DualRand::saveStatus(const char filename[]) const {
134  std::ofstream outFile(filename, std::ios::out);
135  if (!outFile.bad()) {
136  outFile << "Uvec\n";
137  std::vector<unsigned long> v = put();
138  #ifdef TRACE_IO
139  std::cout << "Result of v = put() is:\n";
140  #endif
141  for (unsigned int i=0; i<v.size(); ++i) {
142  outFile << v[i] << "\n";
143  #ifdef TRACE_IO
144  std::cout << v[i] << " ";
145  if (i%6==0) std::cout << "\n";
146  #endif
147  }
148  #ifdef TRACE_IO
149  std::cout << "\n";
150  #endif
151  }
152 #ifdef REMOVED
153  int pr=outFile.precision(20);
154  outFile << theSeed << std::endl;
155  tausworthe.put(outFile);
156  integerCong.put(outFile);
157  outFile << std::endl; // This is superfluous but harmless
158  outFile.precision(pr);
159 #endif
160 }
161 
162 void DualRand::restoreStatus(const char filename[]) {
163  std::ifstream inFile(filename, std::ios::in);
164  if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) {
165  std::cerr << " -- Engine state remains unchanged\n";
166  return;
167  }
168  if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) {
169  std::vector<unsigned long> v;
170  unsigned long xin;
171  for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
172  inFile >> xin;
173  #ifdef TRACE_IO
174  std::cout << "ivec = " << ivec << " xin = " << xin << " ";
175  if (ivec%3 == 0) std::cout << "\n";
176  #endif
177  if (!inFile) {
178  inFile.clear(std::ios::badbit | inFile.rdstate());
179  std::cerr << "\nDualRand state (vector) description improper."
180  << "\nrestoreStatus has failed."
181  << "\nInput stream is probably mispositioned now." << std::endl;
182  return;
183  }
184  v.push_back(xin);
185  }
186  getState(v);
187  return;
188  }
189 
190  if (!inFile.bad()) {
191 // inFile >> theSeed; removed -- encompased by possibleKeywordInput
192  tausworthe.get(inFile);
193  integerCong.get(inFile);
194  }
195 }
196 
197 void DualRand::showStatus() const {
198  int pr=std::cout.precision(20);
199  std::cout << std::endl;
200  std::cout << "-------- DualRand engine status ---------"
201  << std::endl;
202  std::cout << "Initial seed = " << theSeed << std::endl;
203  std::cout << "Tausworthe generator = " << std::endl;
204  tausworthe.put(std::cout);
205  std::cout << "\nIntegerCong generator = " << std::endl;
206  integerCong.put(std::cout);
207  std::cout << std::endl << "-----------------------------------------"
208  << std::endl;
209  std::cout.precision(pr);
210 }
211 
212 DualRand::operator float() {
213  return (float) ( (integerCong ^ tausworthe) * twoToMinus_32()
214  + nearlyTwoToMinus_54() );
215  // add this so that zero never happens
216 }
217 
218 DualRand::operator unsigned int() {
219  return (integerCong ^ tausworthe) & 0xffffffff;
220 }
221 
222 std::ostream & DualRand::put(std::ostream & os) const {
223  char beginMarker[] = "DualRand-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 #ifdef REMOVED
231  char endMarker[] = "DualRand-end";
232  int pr=os.precision(20);
233  os << " " << beginMarker << " ";
234  os << theSeed << " ";
235  tausworthe.put(os);
236  integerCong.put(os);
237  os << " " << endMarker << "\n";
238  os.precision(pr);
239  return os;
240 #endif
241 }
242 
243 std::vector<unsigned long> DualRand::put () const {
244  std::vector<unsigned long> v;
245  v.push_back (engineIDulong<DualRand>());
246  tausworthe.put(v);
247  integerCong.put(v);
248  return v;
249 }
250 
251 std::istream & DualRand::get(std::istream & is) {
252  char beginMarker [MarkerLen];
253  is >> std::ws;
254  is.width(MarkerLen); // causes the next read to the char* to be <=
255  // that many bytes, INCLUDING A TERMINATION \0
256  // (Stroustrup, section 21.3.2)
257  is >> beginMarker;
258  if (strcmp(beginMarker,"DualRand-begin")) {
259  is.clear(std::ios::badbit | is.rdstate());
260  std::cerr << "\nInput mispositioned or"
261  << "\nDualRand state description missing or"
262  << "\nwrong engine type found." << std::endl;
263  return is;
264  }
265  return getState(is);
266 }
267 
268 std::string DualRand::beginTag ( ) {
269  return "DualRand-begin";
270 }
271 
272 std::istream & DualRand::getState ( std::istream & is ) {
273  if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) {
274  std::vector<unsigned long> v;
275  unsigned long uu;
276  for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
277  is >> uu;
278  if (!is) {
279  is.clear(std::ios::badbit | is.rdstate());
280  std::cerr << "\nDualRand state (vector) description improper."
281  << "\ngetState() has failed."
282  << "\nInput stream is probably mispositioned now." << std::endl;
283  return is;
284  }
285  v.push_back(uu);
286  }
287  getState(v);
288  return (is);
289  }
290 
291 // is >> theSeed; Removed, encompassed by possibleKeywordInput()
292 
293  char endMarker [MarkerLen];
294  tausworthe.get(is);
295  integerCong.get(is);
296  is >> std::ws;
297  is.width(MarkerLen);
298  is >> endMarker;
299  if (strcmp(endMarker,"DualRand-end")) {
300  is.clear(std::ios::badbit | is.rdstate());
301  std::cerr << "DualRand state description incomplete."
302  << "\nInput stream is probably mispositioned now." << std::endl;
303  return is;
304  }
305  return is;
306 }
307 
308 bool DualRand::get(const std::vector<unsigned long> & v) {
309  if ((v[0] & 0xffffffffUL) != engineIDulong<DualRand>()) {
310  std::cerr <<
311  "\nDualRand get:state vector has wrong ID word - state unchanged\n";
312  return false;
313  }
314  if (v.size() != VECTOR_STATE_SIZE) {
315  std::cerr << "\nDualRand get:state vector has wrong size: "
316  << v.size() << " - state unchanged\n";
317  return false;
318  }
319  return getState(v);
320 }
321 
322 bool DualRand::getState (const std::vector<unsigned long> & v) {
323  std::vector<unsigned long>::const_iterator iv = v.begin()+1;
324  if (!tausworthe.get(iv)) return false;
325  if (!integerCong.get(iv)) return false;
326  if (iv != v.end()) {
327  std::cerr <<
328  "\nDualRand get:state vector has wrong size: " << v.size()
329  << "\n Apparently " << iv-v.begin() << " words were consumed\n";
330  return false;
331  }
332  return true;
333 }
334 
335 DualRand::Tausworthe::Tausworthe() {
336  words[0] = 1234567;
337  for (wordIndex = 1; wordIndex < 4; ++wordIndex) {
338  words[wordIndex] = 69607 * words[wordIndex-1] + 54329;
339  }
340 }
341 
342 DualRand::Tausworthe::Tausworthe(unsigned int seed) {
343  words[0] = seed;
344  for (wordIndex = 1; wordIndex < 4; ++wordIndex) {
345  words[wordIndex] = 69607 * words[wordIndex-1] + 54329;
346  }
347 }
348 
349 DualRand::Tausworthe::operator unsigned int() {
350 
351 // Mathematically: Consider a sequence of bits b[n]. Repeatedly form
352 // b[0]' = b[127] ^ b[97]; b[n]' = b[n-1]. This sequence will have a very
353 // long period (2**127-1 according to Tausworthe's work).
354 
355 // The actual method used relies on the fact that the operations needed to
356 // form bit 0 for up to 96 iterations never depend on the results of the
357 // previous ones. So you can actually compute many bits at once. In fact
358 // you can compute 32 at once -- despite 127 - 97 < 32 -- but 24 was used in
359 // the method used in Canopy, where they wanted only single-precision float
360 // randoms. I will do 32 here.
361 
362 // When you do it this way, this looks disturbingly like the dread lagged XOR
363 // Fibonacci. And indeed, it is a lagged Fibonacii, F(4,3, op) with the op
364 // being the XOR of a combination of shifts of the two numbers. Although
365 // Tausworthe asserted excellent properties, I would be scared to death.
366 // However, the shifting and bit swapping really does randomize this in a
367 // serious way.
368 
369 // Statements have been made to the effect that shift register sequences fail
370 // the parking lot test because they achieve randomness by multiple foldings,
371 // and this produces a characteristic pattern. We observe that in this
372 // specific algorithm, which has a fairly long lever arm, the foldings become
373 // effectively random. This is evidenced by the fact that the generator
374 // does pass the Diehard tests, including the parking lot test.
375 
376 // To avoid shuffling of variables in memory, you either have to use circular
377 // pointers (and those give you ifs, which are also costly) or compute at least
378 // a few iterations at once. We do the latter. Although there is a possible
379 // trade of room for more speed, by computing and saving 256 instead of 128
380 // bits at once, I will stop at this level of optimization.
381 
382 // To remind: Each (32-bit) step takes the XOR of bits [127-96] with bits
383 // [95-64] and places it in bits [0-31]. But in the first step, we designate
384 // word0 as bits [0-31], in the second step, word 1 (since the bits it holds
385 // will no longer be needed), then word 2, then word 3. After this, the
386 // stream contains 128 random bits which we will use as 4 valid 32-bit
387 // random numbers.
388 
389 // Thus at the start of the first step, word[0] contains the newest (used)
390 // 32-bit random, and word[3] the oldest. After four steps, word[0] again
391 // contains the newest (now unused) random, and word[3] the oldest.
392 // Bit 0 of word[0] is logically the newest bit, and bit 31 of word[3]
393 // the oldest.
394 
395  if (wordIndex <= 0) {
396  for (wordIndex = 0; wordIndex < 4; ++wordIndex) {
397  words[wordIndex] = ( (words[(wordIndex+1) & 3] << 1 ) |
398  (words[wordIndex] >> 31) )
399  ^ ( (words[(wordIndex+1) & 3] << 31) |
400  (words[wordIndex] >> 1) );
401  }
402  }
403  return words[--wordIndex] & 0xffffffff;
404 }
405 
406 void DualRand::Tausworthe::put(std::ostream & os) const {
407  char beginMarker[] = "Tausworthe-begin";
408  char endMarker[] = "Tausworthe-end";
409 
410  int pr=os.precision(20);
411  os << " " << beginMarker << " ";
412  for (int i = 0; i < 4; ++i) {
413  os << words[i] << " ";
414  }
415  os << wordIndex;
416  os << " " << endMarker << " ";
417  os << std::endl;
418  os.precision(pr);
419 }
420 
421 void DualRand::Tausworthe::put(std::vector<unsigned long> & v) const {
422  for (int i = 0; i < 4; ++i) {
423  v.push_back(static_cast<unsigned long>(words[i]));
424  }
425  v.push_back(static_cast<unsigned long>(wordIndex));
426 }
427 
428 void DualRand::Tausworthe::get(std::istream & is) {
429  char beginMarker [MarkerLen];
430  char endMarker [MarkerLen];
431 
432  is >> std::ws;
433  is.width(MarkerLen); // causes the next read to the char* to be <=
434  // that many bytes, INCLUDING A TERMINATION \0
435  // (Stroustrup, section 21.3.2)
436  is >> beginMarker;
437  if (strcmp(beginMarker,"Tausworthe-begin")) {
438  is.clear(std::ios::badbit | is.rdstate());
439  std::cerr << "\nInput mispositioned or"
440  << "\nTausworthe state description missing or"
441  << "\nwrong engine type found." << std::endl;
442  }
443  for (int i = 0; i < 4; ++i) {
444  is >> words[i];
445  }
446  is >> wordIndex;
447  is >> std::ws;
448  is.width(MarkerLen);
449  is >> endMarker;
450  if (strcmp(endMarker,"Tausworthe-end")) {
451  is.clear(std::ios::badbit | is.rdstate());
452  std::cerr << "\nTausworthe state description incomplete."
453  << "\nInput stream is probably mispositioned now." << std::endl;
454  }
455 }
456 
457 bool
458 DualRand::Tausworthe::get(std::vector<unsigned long>::const_iterator & iv){
459  for (int i = 0; i < 4; ++i) {
460  words[i] = *iv++;
461  }
462  wordIndex = *iv++;
463  return true;
464 }
465 
466 DualRand::IntegerCong::IntegerCong()
467 : state((unsigned int)3758656018U),
468  multiplier(66565),
469  addend(12341)
470 {
471 }
472 
473 DualRand::IntegerCong::IntegerCong(unsigned int seed, int streamNumber)
474 : state(seed),
475  multiplier(65536 + 1024 + 5 + (8 * 1017 * streamNumber)),
476  addend(12341)
477 {
478  // As to the multiplier, the following comment was made:
479  // We want our multipliers larger than 2^16, and equal to
480  // 1 mod 4 (for max. period), but not equal to 1 mod 8
481  // (for max. potency -- the smaller and higher dimension the
482  // stripes, the better)
483 
484  // All of these will have fairly long periods. Depending on the choice
485  // of stream number, some of these will be quite decent when considered
486  // as independent random engines, while others will be poor. Thus these
487  // should not be used as stand-alone engines; but when combined with a
488  // generator known to be good, they allow creation of millions of good
489  // independent streams, without fear of two streams accidentally hitting
490  // nearby places in the good random sequence.
491 }
492 
493 DualRand::IntegerCong::operator unsigned int() {
494  return state = (state * multiplier + addend) & 0xffffffff;
495 }
496 
497 void DualRand::IntegerCong::put(std::ostream & os) const {
498  char beginMarker[] = "IntegerCong-begin";
499  char endMarker[] = "IntegerCong-end";
500 
501  int pr=os.precision(20);
502  os << " " << beginMarker << " ";
503  os << state << " " << multiplier << " " << addend;
504  os << " " << endMarker << " ";
505  os << std::endl;
506  os.precision(pr);
507 }
508 
509 void DualRand::IntegerCong::put(std::vector<unsigned long> & v) const {
510  v.push_back(static_cast<unsigned long>(state));
511  v.push_back(static_cast<unsigned long>(multiplier));
512  v.push_back(static_cast<unsigned long>(addend));
513 }
514 
515 void DualRand::IntegerCong::get(std::istream & is) {
516  char beginMarker [MarkerLen];
517  char endMarker [MarkerLen];
518 
519  is >> std::ws;
520  is.width(MarkerLen); // causes the next read to the char* to be <=
521  // that many bytes, INCLUDING A TERMINATION \0
522  // (Stroustrup, section 21.3.2)
523  is >> beginMarker;
524  if (strcmp(beginMarker,"IntegerCong-begin")) {
525  is.clear(std::ios::badbit | is.rdstate());
526  std::cerr << "\nInput mispositioned or"
527  << "\nIntegerCong state description missing or"
528  << "\nwrong engine type found." << std::endl;
529  }
530  is >> state >> multiplier >> addend;
531  is >> std::ws;
532  is.width(MarkerLen);
533  is >> endMarker;
534  if (strcmp(endMarker,"IntegerCong-end")) {
535  is.clear(std::ios::badbit | is.rdstate());
536  std::cerr << "\nIntegerCong state description incomplete."
537  << "\nInput stream is probably mispositioned now." << std::endl;
538  }
539 }
540 
541 bool
542 DualRand::IntegerCong::get(std::vector<unsigned long>::const_iterator & iv) {
543  state = *iv++;
544  multiplier = *iv++;
545  addend = *iv++;
546  return true;
547 }
548 
549 } // namespace CLHEP
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
virtual ~DualRand()
Definition: DualRand.cc:105
static const unsigned int VECTOR_STATE_SIZE
void flatArray(const int size, double *vect)
Definition: DualRand.cc:116
static std::string engineName()
void restoreStatus(const char filename[]="DualRand.conf")
Definition: DualRand.cc:162
static double nearlyTwoToMinus_54()
void saveStatus(const char filename[]="DualRand.conf") const
Definition: DualRand.cc:133
static double twoToMinus_32()
virtual std::istream & get(std::istream &is)
Definition: DualRand.cc:251
double flat()
Definition: DualRand.cc:107
static std::string beginTag()
Definition: DualRand.cc:268
virtual std::istream & getState(std::istream &is)
Definition: DualRand.cc:272
std::vector< unsigned long > put() const
Definition: DualRand.cc:243
void setSeeds(const long *seeds, int)
Definition: DualRand.cc:128
void showStatus() const
Definition: DualRand.cc:197
static double twoToMinus_53()
static bool checkFile(std::istream &file, const std::string &filename, const std::string &classname, const std::string &methodname)
Definition: RandomEngine.cc:46
std::string name() const
Definition: DualRand.cc:63
void setSeed(long seed, int)
Definition: DualRand.cc:122