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