Geant4 Cross Reference |
1 // -*- C++ -*- 2 // 3 // ----------------------------------------------------------------------- 4 // HEP Random 5 // --- RanluxEngine --- 6 // class implementation file 7 // ----------------------------------------------------------------------- 8 // This file is part of Geant4 (simulation toolkit for HEP). 9 // 10 // Ranlux random number generator originally implemented in FORTRAN77 11 // by Fred James as part of the MATHLIB HEP library. 12 // 'RanluxEngine' is designed to fit into the CLHEP random number 13 // class structure. 14 15 // =============================================================== 16 // Adeyemi Adesanya - Created: 6th November 1995 17 // Gabriele Cosmo - Adapted & Revised: 22nd November 1995 18 // Adeyemi Adesanya - Added setSeeds() method: 2nd February 1996 19 // Gabriele Cosmo - Added flatArray() method: 8th February 1996 20 // - Minor corrections: 31st October 1996 21 // - Added methods for engine status: 19th November 1996 22 // - Fixed bug in setSeeds(): 15th September 1997 23 // J.Marraffino - Added stream operators and related constructor. 24 // Added automatic seed selection from seed table and 25 // engine counter: 14th Feb 1998 26 // - Fixed bug: setSeeds() requires a zero terminated 27 // array of seeds: 19th Feb 1998 28 // Ken Smith - Added conversion operators: 6th Aug 1998 29 // J. Marraffino - Remove dependence on hepString class 13 May 1999 30 // M. Fischler - In restore, checkFile for file not found 03 Dec 2004 31 // M. Fischler - Methods put, getfor instance save/restore 12/8/04 32 // M. Fischler - split get() into tag validation and 33 // getState() for anonymous restores 12/27/04 34 // M. Fischler - put/get for vectors of ulongs 3/14/05 35 // M. Fischler - State-saving using only ints, for portability 4/12/05 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 selection 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 room to hold a begin or end marker. 62 63 std::string RanluxEngine::name() const {return "RanluxEngine";} 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/maxIndex)); 88 int curIndex = std::abs(int(numEngines%maxIndex)); 89 90 long mask = ((cycle & 0x007fffff) << 8); 91 HepRandom::getTheTableSeeds( seedlist, curIndex ); 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 colIndex, int lux) 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 Multiplicative 133 // Congruential generator using formula constants of L'Ecuyer 134 // as described in "A review of pseudorandom number generators" 135 // (Fred James) published in Computer Physics Communications 60 (1990) 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 need to be 'thrown away' 151 // every 24 numbers is set using luxury level variable. 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 to default luxury level 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_multiple * ecuyer_a) 169 - k_multiple * ecuyer_c ; 170 if(next_seed < 0)next_seed += ecuyer_d; 171 int_seed_table[i] = next_seed % int_modulus; 172 } 173 174 for(i = 0;i != 24;i++) 175 float_seed_table[i] = int_seed_table[i] * mantissa_bit_24(); 176 177 i_lag = 23; 178 j_lag = 9; 179 carry = 0. ; 180 181 if( float_seed_table[23] == 0. ) carry = mantissa_bit_24(); 182 183 count24 = 0; 184 } 185 186 void RanluxEngine::setSeeds(const long *seeds, int lux) { 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 need to be 'thrown away' 211 // every 24 numbers is set using luxury level variable. 212 213 if( (lux > 4)||(lux < 0) ){ 214 if(lux >= 24){ 215 nskip = lux - 24; 216 }else{ 217 nskip = lux_levels[3]; // corresponds to default luxury level 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_modulus; 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_multiple * ecuyer_a) 234 - k_multiple * ecuyer_c ; 235 if(next_seed < 0)next_seed += ecuyer_d; 236 int_seed_table[i] = next_seed % int_modulus; 237 } 238 } 239 240 for(i = 0;i != 24;i++) 241 float_seed_table[i] = int_seed_table[i] * mantissa_bit_24(); 242 243 i_lag = 23; 244 j_lag = 9; 245 carry = 0. ; 246 247 if( float_seed_table[23] == 0. ) carry = mantissa_bit_24(); 248 249 count24 = 0; 250 } 251 252 void RanluxEngine::saveStatus( const char filename[] ) const 253 { 254 std::ofstream outFile( filename, std::ios::out ) ; 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 filename[] ) 265 { 266 std::ifstream inFile( filename, std::ios::in); 267 if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) { 268 std::cerr << " -- Engine state remains unchanged\n"; 269 return; 270 } 271 if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) { 272 std::vector<unsigned long> v; 273 unsigned long xin; 274 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) { 275 inFile >> xin; 276 if (!inFile) { 277 inFile.clear(std::ios::badbit | inFile.rdstate()); 278 std::cerr << "\nRanluxEngine state (vector) description improper." 279 << "\nrestoreStatus has failed." 280 << "\nInput stream is probably mispositioned now." << std::endl; 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 -- encompased by possibleKeywordInput 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 status ---------" << std::endl; 303 std::cout << " Initial seed = " << theSeed << std::endl; 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_lag = " << j_lag << std::endl; 309 std::cout << " carry = " << carry << ", count24 = " << count24 << std::endl; 310 std::cout << " luxury = " << luxury << " nskip = " << nskip << std::endl; 311 std::cout << "----------------------------------------" << std::endl; 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_table[i_lag] - carry; 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_table[j_lag]; 336 if( uni == 0) uni = mantissa_bit_24() * mantissa_bit_24(); 337 } 338 next_random = uni; 339 count24 ++; 340 341 // every 24th number generation, several random numbers are generated 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_seed_table[i_lag] - carry; 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, double* vect) 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_table[i_lag] - carry; 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_table[j_lag]; 388 if( uni == 0) uni = mantissa_bit_24() * mantissa_bit_24(); 389 } 390 next_random = uni; 391 vect[index] = (double)next_random; 392 count24 ++; 393 394 // every 24th number generation, several random numbers are generated 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] - float_seed_table[i_lag] - carry; 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_bit_32()) & 0xffffffff) | 427 (((unsigned int)(float_seed_table[i_lag]*exponent_bit_32())>>16) & 0xff); 428 // needed because Ranlux doesn't fill all bits of the double 429 // which therefore doesn't fill all bits of the integer. 430 } 431 432 std::ostream & RanluxEngine::put ( std::ostream& os ) const 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 () const { 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_table[i]/mantissa_bit_24())); 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/mantissa_bit_24())); 453 v.push_back(static_cast<unsigned long>(count24)); 454 v.push_back(static_cast<unsigned long>(luxury)); 455 v.push_back(static_cast<unsigned long>(nskip)); 456 return v; 457 } 458 459 std::istream & RanluxEngine::get ( std::istream& is ) 460 { 461 char beginMarker [MarkerLen]; 462 is >> std::ws; 463 is.width(MarkerLen); // causes the next read to the char* to be <= 464 // that many bytes, INCLUDING A TERMINATION \0 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 mispositioned or" 470 << "\nRanluxEngine state description missing or" 471 << "\nwrong engine type found." << std::endl; 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::istream& is ) 482 { 483 if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) { 484 std::vector<unsigned long> v; 485 unsigned long uu; 486 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) { 487 is >> uu; 488 if (!is) { 489 is.clear(std::ios::badbit | is.rdstate()); 490 std::cerr << "\nRanluxEngine state (vector) description improper." 491 << "\ngetState() has failed." 492 << "\nInput stream is probably mispositioned now." << std::endl; 493 return is; 494 } 495 v.push_back(uu); 496 } 497 getState(v); 498 return (is); 499 } 500 501 // is >> theSeed; Removed, encompassed by possibleKeywordInput() 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 description incomplete." 516 << "\nInput stream is probably mispositioned now." << std::endl; 517 return is; 518 } 519 return is; 520 } 521 522 bool RanluxEngine::get (const std::vector<unsigned long> & v) { 523 if ((v[0] & 0xffffffffUL) != engineIDulong<RanluxEngine>()) { 524 std::cerr << 525 "\nRanluxEngine get:state vector has wrong ID word - state unchanged\n"; 526 return false; 527 } 528 return getState(v); 529 } 530 531 bool RanluxEngine::getState (const std::vector<unsigned long> & v) { 532 if (v.size() != VECTOR_STATE_SIZE ) { 533 std::cerr << 534 "\nRanluxEngine get:state vector has wrong length - state unchanged\n"; 535 return false; 536 } 537 for (int i=0; i<24; ++i) { 538 float_seed_table[i] = v[i+1]*mantissa_bit_24(); 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