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