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