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