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