| Geant4 Cross Reference |
1 // -*- C++ -*- 1
2 //
3 // -------------------------------------------
4 // Hep Random
5 // --- DualRand ---
6 // class implementation file
7 // -------------------------------------------
8 // Exclusive or of a feedback shift register
9 // random number generator. The feedback shi
10 // 127 and 97. The integer congruence genera
11 // multiplier for each stream. The multiplie
12 // full period and maximum "potency" for modu
13 // the combined random number generator is 2^
14 // sequences are different for each stream (n
15 // different place).
16 //
17 // In the original generator used on ACPMAPS:
18 // The feedback shift register generates 24 b
19 // the high order 24 bits of the integer cong
20 // used.
21 //
22 // Instead, to conform with more modern engin
23 // 32 bits at a time and use the full 32 bits
24 // generator.
25 //
26 // References:
27 // Knuth
28 // Tausworthe
29 // Golomb
30 //============================================
31 // Ken Smith - Removed pow() from flat()
32 // - Added conversion operators
33 // J. Marraffino - Added some explicit casts
34 // machines where sizeof(int)
35 // M. Fischler - Modified constructors taki
36 // depend on numEngines (same seeds sho
37 // produce same sequences). Default st
38 // depends on numEngines. 14 Sep 1
39 // - Modified use of the variou
40 // to avoid per-instance spac
41 // correct the rounding proce
42 // J. Marraffino - Remove dependence on hepSt
43 // M. Fischler - Put endl at end of a save
44 // M. Fischler - In restore, checkFile for
45 // M. Fischler - methods for distrib. insta
46 // M. Fischler - split get() into tag valid
47 // getState() for anonymous r
48 // Mark Fischler - methods for vector save/re
49 // M. Fischler - State-saving using only in
50 //
51 //============================================
52
53 #include "CLHEP/Random/DualRand.h"
54 #include "CLHEP/Random/engineIDulong.h"
55 #include "CLHEP/Utility/atomic_int.h"
56
57 #include <atomic>
58 #include <ostream>
59 #include <string.h> // for strcmp
60 #include <vector>
61 #include <iostream>
62
63 namespace CLHEP {
64
65 namespace {
66 // Number of instances with automatic seed s
67 CLHEP_ATOMIC_INT_TYPE numberOfEngines(0);
68 }
69
70 static const int MarkerLen = 64; // Enough roo
71
72 std::string DualRand::name() const {return "Du
73
74 // The following constructors (excluding the i
75 // the bits of the tausworthe and the starting
76 // congruence based on the seed. In addition,
77 // for the integer congruence based on the str
78 // absolutely independent streams.
79
80 DualRand::DualRand()
81 : HepRandomEngine(),
82 numEngines(numberOfEngines++),
83 tausworthe (1234567 + numEngines + 175321),
84 integerCong(69607 * tausworthe + 54329, numE
85 {
86 theSeed = 1234567;
87 }
88
89 DualRand::DualRand(long seed)
90 : HepRandomEngine(),
91 numEngines(0),
92 tausworthe ((unsigned int)seed + 175321),
93 integerCong(69607 * tausworthe + 54329, 8043
94 {
95 theSeed = seed;
96 }
97
98 DualRand::DualRand(std::istream & is)
99 : HepRandomEngine(),
100 numEngines(0)
101 {
102 is >> *this;
103 }
104
105 DualRand::DualRand(int rowIndex, int colIndex)
106 : HepRandomEngine(),
107 numEngines(0),
108 tausworthe (rowIndex + 1000 * colIndex + 853
109 integerCong(69607 * tausworthe + 54329, 1123
110 {
111 theSeed = rowIndex;
112 }
113
114 DualRand::~DualRand() { }
115
116 double DualRand::flat() {
117 unsigned int ic ( integerCong );
118 unsigned int t ( tausworthe );
119 return ( (t ^ ic) * twoToMinus_32() +
120 (t >> 11) * twoToMinus_53() +
121 nearlyTwoToMinus_54() // m
122 );
123 }
124
125 void DualRand::flatArray(const int size, doubl
126 for (int i = 0; i < size; ++i) {
127 vect[i] = flat();
128 }
129 }
130
131 void DualRand::setSeed(long seed, int) {
132 theSeed = seed;
133 tausworthe = Tausworthe((unsigned int)seed
134 integerCong = IntegerCong(69607 * tausworthe
135 }
136
137 void DualRand::setSeeds(const long * seeds, in
138 setSeed(seeds ? *seeds : 1234567, 0);
139 theSeeds = seeds;
140 }
141
142 void DualRand::saveStatus(const char filename[
143 std::ofstream outFile(filename, std::ios::ou
144 if (!outFile.bad()) {
145 outFile << "Uvec\n";
146 std::vector<unsigned long> v = put();
147 for (unsigned int i=0; i<v.size(); ++i) {
148 outFile << v[i] << "\n";
149 }
150 }
151 }
152
153 void DualRand::restoreStatus(const char filena
154 std::ifstream inFile(filename, std::ios::in)
155 if (!checkFile ( inFile, filename, engineNam
156 std::cerr << " -- Engine state remains un
157 return;
158 }
159 if ( possibleKeywordInput ( inFile, "Uvec",
160 std::vector<unsigned long> v;
161 unsigned long xin;
162 for (unsigned int ivec=0; ivec < VECTOR_ST
163 inFile >> xin;
164 if (!inFile) {
165 inFile.clear(std::ios::badbit | inFile
166 std::cerr << "\nDualRand state (vector
167 << "\nrestoreStatus has failed."
168 << "\nInput stream is probably mispos
169 return;
170 }
171 v.push_back(xin);
172 }
173 getState(v);
174 return;
175 }
176
177 if (!inFile.bad()) {
178 // inFile >> theSeed; removed -- encompas
179 tausworthe.get(inFile);
180 integerCong.get(inFile);
181 }
182 }
183
184 void DualRand::showStatus() const {
185 long pr=std::cout.precision(20);
186 std::cout << std::endl;
187 std::cout << "-------- DualRand engi
188 << std::endl;
189 std::cout << "Initial seed = " << t
190 std::cout << "Tausworthe generator = " << s
191 tausworthe.put(std::cout);
192 std::cout << "\nIntegerCong generator = " <<
193 integerCong.put(std::cout);
194 std::cout << std::endl << "-----------------
195 << std::endl;
196 std::cout.precision(pr);
197 }
198
199 DualRand::operator double() {
200 return flat();
201 }
202
203 DualRand::operator float() {
204 return (float) ( (integerCong ^ tausworthe)
205 + nearlyTwoToMinus_54() );
206 // add this so that zero never happe
207 }
208
209 DualRand::operator unsigned int() {
210 return (integerCong ^ tausworthe) & 0xffffff
211 }
212
213 std::ostream & DualRand::put(std::ostream & os
214 char beginMarker[] = "DualRand-begin";
215 os << beginMarker << "\nUvec\n";
216 std::vector<unsigned long> v = put();
217 for (unsigned int i=0; i<v.size(); ++i) {
218 os << v[i] << "\n";
219 }
220 return os;
221 }
222
223 std::vector<unsigned long> DualRand::put () co
224 std::vector<unsigned long> v;
225 v.push_back (engineIDulong<DualRand>());
226 tausworthe.put(v);
227 integerCong.put(v);
228 return v;
229 }
230
231 std::istream & DualRand::get(std::istream & is
232 char beginMarker [MarkerLen];
233 is >> std::ws;
234 is.width(MarkerLen); // causes the next rea
235 // that many bytes, INCLUDING A TERMINAT
236 // (Stroustrup, section 21.3.2)
237 is >> beginMarker;
238 if (strcmp(beginMarker,"DualRand-begin")) {
239 is.clear(std::ios::badbit | is.rdstate());
240 std::cerr << "\nInput mispositioned or"
241 << "\nDualRand state description missi
242 << "\nwrong engine type found." << std
243 return is;
244 }
245 return getState(is);
246 }
247
248 std::string DualRand::beginTag ( ) {
249 return "DualRand-begin";
250 }
251
252 std::istream & DualRand::getState ( std::istre
253 if ( possibleKeywordInput ( is, "Uvec", theS
254 std::vector<unsigned long> v;
255 unsigned long uu;
256 for (unsigned int ivec=0; ivec < VECTOR_ST
257 is >> uu;
258 if (!is) {
259 is.clear(std::ios::badbit | is.rdstate
260 std::cerr << "\nDualRand state (vector
261 << "\ngetState() has failed."
262 << "\nInput stream is probably mispos
263 return is;
264 }
265 v.push_back(uu);
266 }
267 getState(v);
268 return (is);
269 }
270
271 // is >> theSeed; Removed, encompassed by po
272
273 char endMarker [MarkerLen];
274 tausworthe.get(is);
275 integerCong.get(is);
276 is >> std::ws;
277 is.width(MarkerLen);
278 is >> endMarker;
279 if (strcmp(endMarker,"DualRand-end")) {
280 is.clear(std::ios::badbit | is.rdstate());
281 std::cerr << "DualRand state description i
282 << "\nInput stream is probably misposi
283 return is;
284 }
285 return is;
286 }
287
288 bool DualRand::get(const std::vector<unsigned
289 if ((v[0] & 0xffffffffUL) != engineIDulong<D
290 std::cerr <<
291 "\nDualRand get:state vector has wrong I
292 return false;
293 }
294 if (v.size() != VECTOR_STATE_SIZE) {
295 std::cerr << "\nDualRand get:state vector
296 << v.size() << " - state unchanged\n";
297 return false;
298 }
299 return getState(v);
300 }
301
302 bool DualRand::getState (const std::vector<uns
303 std::vector<unsigned long>::const_iterator i
304 if (!tausworthe.get(iv)) return false;
305 if (!integerCong.get(iv)) return false;
306 if (iv != v.end()) {
307 std::cerr <<
308 "\nDualRand get:state vector has wrong s
309 << "\n Apparently " << iv-v.begin()
310 return false;
311 }
312 return true;
313 }
314
315 DualRand::Tausworthe::Tausworthe() {
316 words[0] = 1234567;
317 for (wordIndex = 1; wordIndex < 4; ++wordInd
318 words[wordIndex] = 69607 * words[wordIndex
319 }
320 }
321
322 DualRand::Tausworthe::Tausworthe(unsigned int
323 words[0] = seed;
324 for (wordIndex = 1; wordIndex < 4; ++wordInd
325 words[wordIndex] = 69607 * words[wordIndex
326 }
327 }
328
329 DualRand::Tausworthe::operator unsigned int()
330
331 // Mathematically: Consider a sequence of bit
332 // b[0]' = b[127] ^ b[97]; b[n]' = b[n-1]. Th
333 // long period (2**127-1 according to Tauswort
334
335 // The actual method used relies on the fact t
336 // form bit 0 for up to 96 iterations never de
337 // previous ones. So you can actually compute
338 // you can compute 32 at once -- despite 127 -
339 // the method used in Canopy, where they wante
340 // randoms. I will do 32 here.
341
342 // When you do it this way, this looks disturb
343 // Fibonacci. And indeed, it is a lagged Fibo
344 // being the XOR of a combination of shifts of
345 // Tausworthe asserted excellent properties, I
346 // However, the shifting and bit swapping real
347 // serious way.
348
349 // Statements have been made to the effect tha
350 // the parking lot test because they achieve r
351 // and this produces a characteristic pattern.
352 // specific algorithm, which has a fairly long
353 // effectively random. This is evidenced by t
354 // does pass the Diehard tests, including the
355
356 // To avoid shuffling of variables in memory,
357 // pointers (and those give you ifs, which are
358 // a few iterations at once. We do the latter
359 // trade of room for more speed, by computing
360 // bits at once, I will stop at this level of
361
362 // To remind: Each (32-bit) step takes the XO
363 // [95-64] and places it in bits [0-31]. But
364 // word0 as bits [0-31], in the second step, w
365 // will no longer be needed), then word 2, the
366 // stream contains 128 random bits which we wi
367 // random numbers.
368
369 // Thus at the start of the first step, word[0
370 // 32-bit random, and word[3] the oldest. Af
371 // contains the newest (now unused) random, an
372 // Bit 0 of word[0] is logically the newest bi
373 // the oldest.
374
375 if (wordIndex <= 0) {
376 for (wordIndex = 0; wordIndex < 4; ++wordI
377 words[wordIndex] = ( (words[(wordIndex+1
378 (words[word
379 ^ ( (words[(wordIndex+1
380 (words[word
381 }
382 }
383 return words[--wordIndex] & 0xffffffff;
384 }
385
386 void DualRand::Tausworthe::put(std::ostream &
387 char beginMarker[] = "Tausworthe-begin";
388 char endMarker[] = "Tausworthe-end";
389
390 long pr=os.precision(20);
391 os << " " << beginMarker << " ";
392 for (int i = 0; i < 4; ++i) {
393 os << words[i] << " ";
394 }
395 os << wordIndex;
396 os << " " << endMarker << " ";
397 os << std::endl;
398 os.precision(pr);
399 }
400
401 void DualRand::Tausworthe::put(std::vector<uns
402 for (int i = 0; i < 4; ++i) {
403 v.push_back(static_cast<unsigned long>(wor
404 }
405 v.push_back(static_cast<unsigned long>(wordI
406 }
407
408 void DualRand::Tausworthe::get(std::istream &
409 char beginMarker [MarkerLen];
410 char endMarker [MarkerLen];
411
412 is >> std::ws;
413 is.width(MarkerLen); // causes the next rea
414 // that many bytes, INCLUDING A TERMINAT
415 // (Stroustrup, section 21.3.2)
416 is >> beginMarker;
417 if (strcmp(beginMarker,"Tausworthe-begin"))
418 is.clear(std::ios::badbit | is.rdstate());
419 std::cerr << "\nInput mispositioned or"
420 << "\nTausworthe state description mis
421 << "\nwrong engine type found." << std
422 }
423 for (int i = 0; i < 4; ++i) {
424 is >> words[i];
425 }
426 is >> wordIndex;
427 is >> std::ws;
428 is.width(MarkerLen);
429 is >> endMarker;
430 if (strcmp(endMarker,"Tausworthe-end")) {
431 is.clear(std::ios::badbit | is.rdstate());
432 std::cerr << "\nTausworthe state descripti
433 << "\nInput stream is probably misposi
434 }
435 }
436
437 bool
438 DualRand::Tausworthe::get(std::vector<unsigned
439 for (int i = 0; i < 4; ++i) {
440 words[i] = (unsigned int)*iv++;
441 }
442 wordIndex = (int)*iv++;
443 return true;
444 }
445
446 DualRand::IntegerCong::IntegerCong()
447 : state((unsigned int)3758656018U),
448 multiplier(66565),
449 addend(12341)
450 {
451 }
452
453 DualRand::IntegerCong::IntegerCong(unsigned in
454 : state(seed),
455 multiplier(65536 + 1024 + 5 + (8 * 1017 * st
456 addend(12341)
457 {
458 // As to the multiplier, the following comme
459 // We want our multipliers larger than 2^16,
460 // 1 mod 4 (for max. period), but not equal
461 // (for max. potency -- the smaller and high
462 // stripes, the better)
463
464 // All of these will have fairly long period
465 // of stream number, some of these will be q
466 // as independent random engines, while othe
467 // should not be used as stand-alone engines
468 // generator known to be good, they allow cr
469 // independent streams, without fear of two
470 // nearby places in the good random sequence
471 }
472
473 DualRand::IntegerCong::operator unsigned int()
474 return state = (state * multiplier + addend)
475 }
476
477 void DualRand::IntegerCong::put(std::ostream &
478 char beginMarker[] = "IntegerCong-begin";
479 char endMarker[] = "IntegerCong-end";
480
481 long pr=os.precision(20);
482 os << " " << beginMarker << " ";
483 os << state << " " << multiplier << " " << a
484 os << " " << endMarker << " ";
485 os << std::endl;
486 os.precision(pr);
487 }
488
489 void DualRand::IntegerCong::put(std::vector<un
490 v.push_back(static_cast<unsigned long>(state
491 v.push_back(static_cast<unsigned long>(multi
492 v.push_back(static_cast<unsigned long>(adden
493 }
494
495 void DualRand::IntegerCong::get(std::istream &
496 char beginMarker [MarkerLen];
497 char endMarker [MarkerLen];
498
499 is >> std::ws;
500 is.width(MarkerLen); // causes the next rea
501 // that many bytes, INCLUDING A TERMINAT
502 // (Stroustrup, section 21.3.2)
503 is >> beginMarker;
504 if (strcmp(beginMarker,"IntegerCong-begin"))
505 is.clear(std::ios::badbit | is.rdstate());
506 std::cerr << "\nInput mispositioned or"
507 << "\nIntegerCong state description mi
508 << "\nwrong engine type found." << std
509 }
510 is >> state >> multiplier >> addend;
511 is >> std::ws;
512 is.width(MarkerLen);
513 is >> endMarker;
514 if (strcmp(endMarker,"IntegerCong-end")) {
515 is.clear(std::ios::badbit | is.rdstate());
516 std::cerr << "\nIntegerCong state descript
517 << "\nInput stream is probably misposi
518 }
519 }
520
521 bool
522 DualRand::IntegerCong::get(std::vector<unsigne
523 state = (unsigned int)*iv++;
524 multiplier = (unsigned int)*iv++;
525 addend = (unsigned int)*iv++;
526 return true;
527 }
528
529 } // namespace CLHEP
530