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