Geant4 Cross Reference |
1 // -*- C++ -*- 1 2 // 3 // ------------------------------------------- 4 // HEP Random 5 // --- RanluxEngine --- 6 // class implementation f 7 // ------------------------------------------- 8 // This file is part of Geant4 (simulation too 9 // 10 // Ranlux random number generator originally i 11 // by Fred James as part of the MATHLIB HEP li 12 // 'RanluxEngine' is designed to fit into the 13 // class structure. 14 15 // =========================================== 16 // Adeyemi Adesanya - Created: 6th November 19 17 // Gabriele Cosmo - Adapted & Revised: 22nd No 18 // Adeyemi Adesanya - Added setSeeds() method: 19 // Gabriele Cosmo - Added flatArray() method: 20 // - Minor corrections: 31st Oc 21 // - Added methods for engine s 22 // - Fixed bug in setSeeds(): 1 23 // J.Marraffino - Added stream operators and 24 // Added automatic seed selec 25 // engine counter: 14th Feb 1 26 // - Fixed bug: setSeeds() requ 27 // array of seeds: 19th Feb 1 28 // Ken Smith - Added conversion operators 29 // J. Marraffino - Remove dependence on hepSt 30 // M. Fischler - In restore, checkFile for 31 // M. Fischler - Methods put, getfor instan 32 // M. Fischler - split get() into tag valid 33 // getState() for anonymous r 34 // M. Fischler - put/get for vectors of ulo 35 // M. Fischler - State-saving using only in 36 // 37 // =========================================== 38 39 #include "CLHEP/Random/Random.h" 40 #include "CLHEP/Random/RanluxEngine.h" 41 #include "CLHEP/Random/engineIDulong.h" 42 #include "CLHEP/Utility/atomic_int.h" 43 44 #include <atomic> 45 #include <cstdlib> // for std::abs(int) 46 #include <iostream> 47 #include <string.h> // for strcmp 48 #include <string> 49 #include <vector> 50 51 namespace CLHEP { 52 53 namespace { 54 // Number of instances with automatic seed s 55 CLHEP_ATOMIC_INT_TYPE numberOfEngines(0); 56 57 // Maximum index into the seed table 58 const int maxIndex = 215; 59 } 60 61 static const int MarkerLen = 64; // Enough roo 62 63 std::string RanluxEngine::name() const {return 64 65 RanluxEngine::RanluxEngine(long seed, int lux) 66 : HepRandomEngine() 67 { 68 long seedlist[2]={0,0}; 69 70 luxury = lux; 71 setSeed(seed, luxury); 72 73 // setSeeds() wants a zero terminated array 74 seedlist[0]=theSeed; 75 seedlist[1]=0; 76 setSeeds(seedlist, luxury); 77 } 78 79 RanluxEngine::RanluxEngine() 80 : HepRandomEngine() 81 { 82 long seed; 83 long seedlist[2]={0,0}; 84 85 luxury = 3; 86 int numEngines = numberOfEngines++; 87 int cycle = std::abs(int(numEngines/maxInde 88 int curIndex = std::abs(int(numEngines%maxI 89 90 long mask = ((cycle & 0x007fffff) << 8); 91 HepRandom::getTheTableSeeds( seedlist, curI 92 seed = seedlist[0]^mask; 93 setSeed(seed, luxury); 94 95 // setSeeds() wants a zero terminated array 96 seedlist[0]=theSeed; 97 seedlist[1]=0; 98 setSeeds(seedlist, luxury); 99 } 100 101 RanluxEngine::RanluxEngine(int rowIndex, int c 102 : HepRandomEngine() 103 { 104 long seed; 105 long seedlist[2]={0,0}; 106 107 luxury = lux; 108 int cycle = std::abs(int(rowIndex/maxIndex) 109 int row = std::abs(int(rowIndex%maxIndex)); 110 int col = std::abs(int(colIndex%2)); 111 long mask = (( cycle & 0x000007ff ) << 20 ) 112 HepRandom::getTheTableSeeds( seedlist, row 113 seed = ( seedlist[col] )^mask; 114 setSeed(seed, luxury); 115 116 // setSeeds() wants a zero terminated array 117 seedlist[0]=theSeed; 118 seedlist[1]=0; 119 setSeeds(seedlist, luxury); 120 } 121 122 RanluxEngine::RanluxEngine( std::istream& is ) 123 : HepRandomEngine() 124 { 125 is >> *this; 126 } 127 128 RanluxEngine::~RanluxEngine() {} 129 130 void RanluxEngine::setSeed(long seed, int lux) 131 132 // The initialisation is carried out using a M 133 // Congruential generator using formula consta 134 // as described in "A review of pseudorandom n 135 // (Fred James) published in Computer Physics 136 // pages 329-344 137 138 const int ecuyer_a = 53668; 139 const int ecuyer_b = 40014; 140 const int ecuyer_c = 12211; 141 const int ecuyer_d = 2147483563; 142 143 const int lux_levels[5] = {0,24,73,199,365}; 144 145 long int_seed_table[24]; 146 long next_seed = seed; 147 long k_multiple; 148 int i; 149 150 // number of additional random numbers that ne 151 // every 24 numbers is set using luxury level 152 153 theSeed = seed; 154 if( (lux > 4)||(lux < 0) ){ 155 if(lux >= 24){ 156 nskip = lux - 24; 157 }else{ 158 nskip = lux_levels[3]; // corresponds 159 } 160 }else{ 161 luxury = lux; 162 nskip = lux_levels[luxury]; 163 } 164 165 166 for(i = 0;i != 24;i++){ 167 k_multiple = next_seed / ecuyer_a; 168 next_seed = ecuyer_b * (next_seed - k_mul 169 - k_multiple * ecuyer_c ; 170 if(next_seed < 0)next_seed += ecuyer_d; 171 int_seed_table[i] = next_seed % int_modul 172 } 173 174 for(i = 0;i != 24;i++) 175 float_seed_table[i] = int_seed_table[i] * 176 177 i_lag = 23; 178 j_lag = 9; 179 carry = 0. ; 180 181 if( float_seed_table[23] == 0. ) carry = man 182 183 count24 = 0; 184 } 185 186 void RanluxEngine::setSeeds(const long *seeds, 187 188 const int ecuyer_a = 53668; 189 const int ecuyer_b = 40014; 190 const int ecuyer_c = 12211; 191 const int ecuyer_d = 2147483563; 192 193 const int lux_levels[5] = {0,24,73,199,365} 194 int i; 195 long int_seed_table[24]; 196 long k_multiple,next_seed; 197 const long *seedptr; 198 199 theSeeds = seeds; 200 seedptr = seeds; 201 202 if(seeds == 0){ 203 setSeed(theSeed,lux); 204 theSeeds = &theSeed; 205 return; 206 } 207 208 theSeed = *seeds; 209 210 // number of additional random numbers that ne 211 // every 24 numbers is set using luxury level 212 213 if( (lux > 4)||(lux < 0) ){ 214 if(lux >= 24){ 215 nskip = lux - 24; 216 }else{ 217 nskip = lux_levels[3]; // corresponds 218 } 219 }else{ 220 luxury = lux; 221 nskip = lux_levels[luxury]; 222 } 223 224 for( i = 0;(i != 24)&&(*seedptr != 0);i++){ 225 int_seed_table[i] = *seedptr % int_modul 226 seedptr++; 227 } 228 229 if(i != 24){ 230 next_seed = int_seed_table[i-1]; 231 for(;i != 24;i++){ 232 k_multiple = next_seed / ecuyer_a; 233 next_seed = ecuyer_b * (next_seed - k_ 234 - k_multiple * ecuyer_c ; 235 if(next_seed < 0)next_seed += ecuyer_d 236 int_seed_table[i] = next_seed % int_mo 237 } 238 } 239 240 for(i = 0;i != 24;i++) 241 float_seed_table[i] = int_seed_table[i] * 242 243 i_lag = 23; 244 j_lag = 9; 245 carry = 0. ; 246 247 if( float_seed_table[23] == 0. ) carry = man 248 249 count24 = 0; 250 } 251 252 void RanluxEngine::saveStatus( const char file 253 { 254 std::ofstream outFile( filename, std::ios:: 255 if (!outFile.bad()) { 256 outFile << "Uvec\n"; 257 std::vector<unsigned long> v = put(); 258 for (unsigned int i=0; i<v.size(); ++i) { 259 outFile << v[i] << "\n"; 260 } 261 } 262 } 263 264 void RanluxEngine::restoreStatus( const char f 265 { 266 std::ifstream inFile( filename, std::ios::i 267 if (!checkFile ( inFile, filename, engineNa 268 std::cerr << " -- Engine state remains u 269 return; 270 } 271 if ( possibleKeywordInput ( inFile, "Uvec", 272 std::vector<unsigned long> v; 273 unsigned long xin; 274 for (unsigned int ivec=0; ivec < VECTOR_ST 275 inFile >> xin; 276 if (!inFile) { 277 inFile.clear(std::ios::badbit | inFile 278 std::cerr << "\nRanluxEngine state (ve 279 << "\nrestoreStatus has failed." 280 << "\nInput stream is probably mispos 281 return; 282 } 283 v.push_back(xin); 284 } 285 getState(v); 286 return; 287 } 288 289 if (!inFile.bad() && !inFile.eof()) { 290 // inFile >> theSeed; removed -- encompas 291 for (int i=0; i<24; ++i) 292 inFile >> float_seed_table[i]; 293 inFile >> i_lag; inFile >> j_lag; 294 inFile >> carry; inFile >> count24; 295 inFile >> luxury; inFile >> nskip; 296 } 297 } 298 299 void RanluxEngine::showStatus() const 300 { 301 std::cout << std::endl; 302 std::cout << "--------- Ranlux engine statu 303 std::cout << " Initial seed = " << theSeed 304 std::cout << " float_seed_table[] = "; 305 for (int i=0; i<24; ++i) 306 std::cout << float_seed_table[i] << " "; 307 std::cout << std::endl; 308 std::cout << " i_lag = " << i_lag << ", j_l 309 std::cout << " carry = " << carry << ", cou 310 std::cout << " luxury = " << luxury << " ns 311 std::cout << "----------------------------- 312 } 313 314 double RanluxEngine::flat() { 315 316 float next_random; 317 float uni; 318 int i; 319 320 uni = float_seed_table[j_lag] - float_seed_t 321 if(uni < 0. ){ 322 uni += 1.0; 323 carry = mantissa_bit_24(); 324 }else{ 325 carry = 0.; 326 } 327 328 float_seed_table[i_lag] = uni; 329 i_lag --; 330 j_lag --; 331 if(i_lag < 0) i_lag = 23; 332 if(j_lag < 0) j_lag = 23; 333 334 if( uni < mantissa_bit_12() ){ 335 uni += mantissa_bit_24() * float_seed_tab 336 if( uni == 0) uni = mantissa_bit_24() * m 337 } 338 next_random = uni; 339 count24 ++; 340 341 // every 24th number generation, several rando 342 // and wasted depending upon the luxury level. 343 344 if(count24 == 24 ){ 345 count24 = 0; 346 for( i = 0; i != nskip ; i++){ 347 uni = float_seed_table[j_lag] - float 348 if(uni < 0. ){ 349 uni += 1.0; 350 carry = mantissa_bit_24(); 351 }else{ 352 carry = 0.; 353 } 354 float_seed_table[i_lag] = uni; 355 i_lag --; 356 j_lag --; 357 if(i_lag < 0)i_lag = 23; 358 if(j_lag < 0) j_lag = 23; 359 } 360 } 361 return (double) next_random; 362 } 363 364 void RanluxEngine::flatArray(const int size, d 365 { 366 float next_random; 367 float uni; 368 int i; 369 int index; 370 371 for (index=0; index<size; ++index) { 372 uni = float_seed_table[j_lag] - float_seed 373 if(uni < 0. ){ 374 uni += 1.0; 375 carry = mantissa_bit_24(); 376 }else{ 377 carry = 0.; 378 } 379 380 float_seed_table[i_lag] = uni; 381 i_lag --; 382 j_lag --; 383 if(i_lag < 0) i_lag = 23; 384 if(j_lag < 0) j_lag = 23; 385 386 if( uni < mantissa_bit_12() ){ 387 uni += mantissa_bit_24() * float_seed_t 388 if( uni == 0) uni = mantissa_bit_24() * 389 } 390 next_random = uni; 391 vect[index] = (double)next_random; 392 count24 ++; 393 394 // every 24th number generation, several rando 395 // and wasted depending upon the luxury level. 396 397 if(count24 == 24 ){ 398 count24 = 0; 399 for( i = 0; i != nskip ; i++){ 400 uni = float_seed_table[j_lag] - flo 401 if(uni < 0. ){ 402 uni += 1.0; 403 carry = mantissa_bit_24(); 404 }else{ 405 carry = 0.; 406 } 407 float_seed_table[i_lag] = uni; 408 i_lag --; 409 j_lag --; 410 if(i_lag < 0)i_lag = 23; 411 if(j_lag < 0) j_lag = 23; 412 } 413 } 414 } 415 } 416 417 RanluxEngine::operator double() { 418 return flat(); 419 } 420 421 RanluxEngine::operator float() { 422 return float( flat() ); 423 } 424 425 RanluxEngine::operator unsigned int() { 426 return ((unsigned int)(flat() * exponent_bi 427 (((unsigned int)(float_seed_table[i_l 428 // needed because Ranlux doesn't fill all b 429 // which therefore doesn't fill all bits of 430 } 431 432 std::ostream & RanluxEngine::put ( std::ostrea 433 { 434 char beginMarker[] = "RanluxEngine-begin"; 435 os << beginMarker << "\nUvec\n"; 436 std::vector<unsigned long> v = put(); 437 for (unsigned int i=0; i<v.size(); ++i) { 438 os << v[i] << "\n"; 439 } 440 return os; 441 } 442 443 std::vector<unsigned long> RanluxEngine::put ( 444 std::vector<unsigned long> v; 445 v.push_back (engineIDulong<RanluxEngine>()); 446 for (int i=0; i<24; ++i) { 447 v.push_back 448 (static_cast<unsigned long>(float_seed_t 449 } 450 v.push_back(static_cast<unsigned long>(i_lag 451 v.push_back(static_cast<unsigned long>(j_lag 452 v.push_back(static_cast<unsigned long>(carry 453 v.push_back(static_cast<unsigned long>(count 454 v.push_back(static_cast<unsigned long>(luxur 455 v.push_back(static_cast<unsigned long>(nskip 456 return v; 457 } 458 459 std::istream & RanluxEngine::get ( std::istrea 460 { 461 char beginMarker [MarkerLen]; 462 is >> std::ws; 463 is.width(MarkerLen); // causes the next rea 464 // that many bytes, INCLUDING A TERMINAT 465 // (Stroustrup, section 21.3.2) 466 is >> beginMarker; 467 if (strcmp(beginMarker,"RanluxEngine-begin") 468 is.clear(std::ios::badbit | is.rdstate()) 469 std::cerr << "\nInput stream mispositione 470 << "\nRanluxEngine state description 471 << "\nwrong engine type found." << st 472 return is; 473 } 474 return getState(is); 475 } 476 477 std::string RanluxEngine::beginTag ( ) { 478 return "RanluxEngine-begin"; 479 } 480 481 std::istream & RanluxEngine::getState ( std::i 482 { 483 if ( possibleKeywordInput ( is, "Uvec", theS 484 std::vector<unsigned long> v; 485 unsigned long uu; 486 for (unsigned int ivec=0; ivec < VECTOR_ST 487 is >> uu; 488 if (!is) { 489 is.clear(std::ios::badbit | is.rdstate 490 std::cerr << "\nRanluxEngine state (ve 491 << "\ngetState() has failed." 492 << "\nInput stream is probably mispos 493 return is; 494 } 495 v.push_back(uu); 496 } 497 getState(v); 498 return (is); 499 } 500 501 // is >> theSeed; Removed, encompassed by po 502 503 char endMarker [MarkerLen]; 504 for (int i=0; i<24; ++i) { 505 is >> float_seed_table[i]; 506 } 507 is >> i_lag; is >> j_lag; 508 is >> carry; is >> count24; 509 is >> luxury; is >> nskip; 510 is >> std::ws; 511 is.width(MarkerLen); 512 is >> endMarker; 513 if (strcmp(endMarker,"RanluxEngine-end")) { 514 is.clear(std::ios::badbit | is.rdstate()) 515 std::cerr << "\nRanluxEngine state descri 516 << "\nInput stream is probably mispos 517 return is; 518 } 519 return is; 520 } 521 522 bool RanluxEngine::get (const std::vector<unsi 523 if ((v[0] & 0xffffffffUL) != engineIDulong<R 524 std::cerr << 525 "\nRanluxEngine get:state vector has wro 526 return false; 527 } 528 return getState(v); 529 } 530 531 bool RanluxEngine::getState (const std::vector 532 if (v.size() != VECTOR_STATE_SIZE ) { 533 std::cerr << 534 "\nRanluxEngine get:state vector has wro 535 return false; 536 } 537 for (int i=0; i<24; ++i) { 538 float_seed_table[i] = v[i+1]*mantissa_bit_ 539 } 540 i_lag = (int)v[25]; 541 j_lag = (int)v[26]; 542 carry = v[27]*mantissa_bit_24(); 543 count24 = (int)v[28]; 544 luxury = (int)v[29]; 545 nskip = (int)v[30]; 546 return true; 547 } 548 549 } // namespace CLHEP 550