Geant4 Cross Reference |
>> 1 // $Id:$ 1 // -*- C++ -*- 2 // -*- C++ -*- 2 // 3 // 3 // ------------------------------------------- 4 // ----------------------------------------------------------------------- 4 // HEP Random 5 // HEP Random 5 // --- Ranlux64Engine -- 6 // --- Ranlux64Engine --- 6 // class implementation f 7 // class implementation file 7 // ------------------------------------------- 8 // ----------------------------------------------------------------------- 8 // A double-precision implementation of the Ra 9 // A double-precision implementation of the RanluxEngine generator as 9 // decsribed by the notes of the original ranl 10 // decsribed by the notes of the original ranlux author (Martin Luscher) 10 // 11 // 11 // See the note by Martin Luscher, December 19 12 // See the note by Martin Luscher, December 1997, entitiled 12 // Double-precision implementation of the rand 13 // Double-precision implementation of the random number generator ranlux 13 // 14 // 14 // =========================================== 15 // ======================================================================= 15 // Ken Smith - Initial draft: 14th Jul 19 16 // Ken Smith - Initial draft: 14th Jul 1998 16 // - Removed pow() from flat me 17 // - Removed pow() from flat method 14th Jul 1998 17 // - Added conversion operators 18 // - Added conversion operators: 6th Aug 1998 18 // 19 // 19 // Mark Fischler The following were modified 20 // Mark Fischler The following were modified mostly to make the routine 20 // exactly match the Luscher algorithm in 21 // exactly match the Luscher algorithm in generating 48-bit 21 // randoms: 22 // randoms: 22 // 9/9/98 - Substantial changes in what used 23 // 9/9/98 - Substantial changes in what used to be flat() to match 23 // algorithm in Luscher's ranlxd.c 24 // algorithm in Luscher's ranlxd.c 24 // - Added update() method for 12 numbers 25 // - Added update() method for 12 numbers, making flat() trivial 25 // - Added advance() method to hold the u 26 // - Added advance() method to hold the unrolled loop for update 26 // - Distinction between three forms of s 27 // - Distinction between three forms of seeding such that it 27 // is impossible to get same sequence f 28 // is impossible to get same sequence from different forms - 28 // done by discarding some fraction of 29 // done by discarding some fraction of one macro cycle which 29 // is different for the three cases 30 // is different for the three cases 30 // - Change the misnomer "seed_table" to 31 // - Change the misnomer "seed_table" to the more accurate 31 // "randoms" 32 // "randoms" 32 // - Removed the no longer needed count12 33 // - Removed the no longer needed count12, i_lag, j_lag, etc. 33 // - Corrected seed procedure which had b 34 // - Corrected seed procedure which had been filling bits past 34 // 2^-48. This actually was very bad, 35 // 2^-48. This actually was very bad, invalidating the 35 // number theory behind the proof that 36 // number theory behind the proof that ranlxd is good. 36 // - Addition of 2**(-49) to generated nu 37 // - Addition of 2**(-49) to generated number to prevent zero 37 // from being returned; this does not a 38 // from being returned; this does not affect the sequence 38 // itself. 39 // itself. 39 // - Corrected ecu seeding, which had bee 40 // - Corrected ecu seeding, which had been supplying only 40 // numbers less than 1/2. This is prob 41 // numbers less than 1/2. This is probably moot. 41 // 9/15/98 - Modified use of the various ex 42 // 9/15/98 - Modified use of the various exponents of 2 42 // to avoid per-instance spac 43 // to avoid per-instance space overhead. Note that these 43 // are initialized in setSeed, which EV 44 // are initialized in setSeed, which EVERY constructor 44 // must invoke. 45 // must invoke. 45 // J. Marraffino - Remove dependence on hepSt 46 // J. Marraffino - Remove dependence on hepString class 13 May 1999 46 // M. Fischler - In restore, checkFile for 47 // M. Fischler - In restore, checkFile for file not found 03 Dec 2004 47 // M. Fischler - put get Methods for distri 48 // M. Fischler - put get Methods for distrib instance save/restore 12/8/04 48 // M. Fischler - split get() into tag valid 49 // M. Fischler - split get() into tag validation and 49 // getState() for anonymous r 50 // getState() for anonymous restores 12/27/04 50 // M. Fischler - put/get for vectors of ulo 51 // M. Fischler - put/get for vectors of ulongs 3/14/05 51 // M. Fischler - State-saving using only in 52 // M. Fischler - State-saving using only ints, for portability 4/12/05 52 // 53 // 53 // =========================================== 54 // ======================================================================= 54 55 55 #include "CLHEP/Random/Random.h" 56 #include "CLHEP/Random/Random.h" 56 #include "CLHEP/Random/Ranlux64Engine.h" 57 #include "CLHEP/Random/Ranlux64Engine.h" 57 #include "CLHEP/Random/engineIDulong.h" 58 #include "CLHEP/Random/engineIDulong.h" 58 #include "CLHEP/Random/DoubConv.h" 59 #include "CLHEP/Random/DoubConv.h" 59 #include "CLHEP/Utility/atomic_int.h" << 60 #include <string.h> // for strcmp 60 << 61 #include <atomic> << 62 #include <cstdlib> // for std::abs(int) 61 #include <cstdlib> // for std::abs(int) 63 #include <iostream> << 64 #include <limits> // for numeric_limits 62 #include <limits> // for numeric_limits 65 #include <string.h> // for strcmp << 66 #include <vector> << 67 63 68 namespace CLHEP { 64 namespace CLHEP { 69 65 70 namespace { << 66 static const int MarkerLen = 64; // Enough room to hold a begin or end marker. 71 // Number of instances with automatic seed s << 72 CLHEP_ATOMIC_INT_TYPE numberOfEngines(0); << 73 67 74 // Maximum index into the seed table << 75 const int maxIndex = 215; << 76 } << 77 68 78 static const int MarkerLen = 64; // Enough roo << 69 // Number of instances with automatic seed selection >> 70 int Ranlux64Engine::numEngines = 0; 79 71 >> 72 // Maximum index into the seed table >> 73 int Ranlux64Engine::maxIndex = 215; 80 74 81 #ifndef WIN32 75 #ifndef WIN32 82 namespace detail { 76 namespace detail { 83 77 84 template< std::size_t n, 78 template< std::size_t n, 85 bool = n < std::size_t(std::numeric_ 79 bool = n < std::size_t(std::numeric_limits<unsigned long>::digits) > 86 struct do_right_shift; 80 struct do_right_shift; 87 template< std::size_t n > 81 template< std::size_t n > 88 struct do_right_shift<n,true> 82 struct do_right_shift<n,true> 89 { 83 { 90 unsigned long operator()(unsigned long value 84 unsigned long operator()(unsigned long value) { return value >> n; } 91 }; 85 }; 92 template< std::size_t n > 86 template< std::size_t n > 93 struct do_right_shift<n,false> 87 struct do_right_shift<n,false> 94 { 88 { 95 unsigned long operator()(unsigned long) { re 89 unsigned long operator()(unsigned long) { return 0ul; } 96 }; 90 }; 97 91 98 template< std::size_t nbits > 92 template< std::size_t nbits > 99 unsigned long rshift( unsigned long value ) 93 unsigned long rshift( unsigned long value ) 100 { return do_right_shift<nbits>()(value); } 94 { return do_right_shift<nbits>()(value); } 101 95 102 } // namespace detail 96 } // namespace detail 103 #endif 97 #endif 104 98 105 std::string Ranlux64Engine::name() const {retu 99 std::string Ranlux64Engine::name() const {return "Ranlux64Engine";} 106 100 107 Ranlux64Engine::Ranlux64Engine() 101 Ranlux64Engine::Ranlux64Engine() 108 : HepRandomEngine() 102 : HepRandomEngine() 109 { 103 { 110 luxury = 1; 104 luxury = 1; 111 int numEngines = numberOfEngines++; << 112 int cycle = std::abs(int(numEngines/maxI 105 int cycle = std::abs(int(numEngines/maxIndex)); 113 int curIndex = std::abs(int(numEngines%maxI 106 int curIndex = std::abs(int(numEngines%maxIndex)); 114 << 107 numEngines +=1; 115 long mask = ((cycle & 0x007fffff) << 8); 108 long mask = ((cycle & 0x007fffff) << 8); 116 long seedlist[2]; 109 long seedlist[2]; 117 HepRandom::getTheTableSeeds( seedlist, curI 110 HepRandom::getTheTableSeeds( seedlist, curIndex ); 118 seedlist[0] ^= mask; 111 seedlist[0] ^= mask; 119 seedlist[1] = 0; 112 seedlist[1] = 0; 120 113 121 setSeeds(seedlist, luxury); 114 setSeeds(seedlist, luxury); 122 advance ( 8 ); // Discard some iteratio 115 advance ( 8 ); // Discard some iterations and ensure that 123 // this sequence won't match one where 116 // this sequence won't match one where seeds 124 // were provided. 117 // were provided. 125 } 118 } 126 119 127 Ranlux64Engine::Ranlux64Engine(long seed, int 120 Ranlux64Engine::Ranlux64Engine(long seed, int lux) 128 : HepRandomEngine() 121 : HepRandomEngine() 129 { 122 { 130 luxury = lux; 123 luxury = lux; 131 long seedlist[2]={seed,0}; 124 long seedlist[2]={seed,0}; 132 setSeeds(seedlist, lux); 125 setSeeds(seedlist, lux); 133 advance ( 2*lux + 1 ); // Discard some it 126 advance ( 2*lux + 1 ); // Discard some iterations to use a different 134 // point in the sequence. 127 // point in the sequence. 135 } 128 } 136 129 137 Ranlux64Engine::Ranlux64Engine(int rowIndex, i 130 Ranlux64Engine::Ranlux64Engine(int rowIndex, int, int lux) 138 : HepRandomEngine() 131 : HepRandomEngine() 139 { 132 { 140 luxury = lux; 133 luxury = lux; 141 int cycle = std::abs(int(rowIndex/maxIndex) 134 int cycle = std::abs(int(rowIndex/maxIndex)); 142 int row = std::abs(int(rowIndex%maxIndex) 135 int row = std::abs(int(rowIndex%maxIndex)); 143 long mask = (( cycle & 0x000007ff ) << 20 ) 136 long mask = (( cycle & 0x000007ff ) << 20 ); 144 long seedlist[2]; 137 long seedlist[2]; 145 HepRandom::getTheTableSeeds( seedlist, row 138 HepRandom::getTheTableSeeds( seedlist, row ); 146 seedlist[0] ^= mask; 139 seedlist[0] ^= mask; 147 seedlist[1]= 0; 140 seedlist[1]= 0; 148 setSeeds(seedlist, lux); 141 setSeeds(seedlist, lux); 149 } 142 } 150 143 151 Ranlux64Engine::Ranlux64Engine( std::istream& 144 Ranlux64Engine::Ranlux64Engine( std::istream& is ) 152 : HepRandomEngine() 145 : HepRandomEngine() 153 { 146 { 154 is >> *this; 147 is >> *this; 155 } 148 } 156 149 157 Ranlux64Engine::~Ranlux64Engine() {} 150 Ranlux64Engine::~Ranlux64Engine() {} 158 151 159 double Ranlux64Engine::flat() { 152 double Ranlux64Engine::flat() { 160 // Luscher improves the speed by computing s 153 // Luscher improves the speed by computing several numbers in a shot, 161 // in a manner similar to that of the Tauswo 154 // in a manner similar to that of the Tausworth in DualRand or the Hurd 162 // engines. Thus, the real work is done in 155 // engines. Thus, the real work is done in update(). Here we merely ensure 163 // that zero, which the algorithm can produc 156 // that zero, which the algorithm can produce, is never returned by flat(). 164 157 165 if (index <= 0) update(); 158 if (index <= 0) update(); 166 return randoms[--index] + twoToMinus_49(); 159 return randoms[--index] + twoToMinus_49(); 167 } 160 } 168 161 169 void Ranlux64Engine::update() { 162 void Ranlux64Engine::update() { 170 // Update the stash of twelve random numbers 163 // Update the stash of twelve random numbers. 171 // When this routione is entered, index is a 164 // When this routione is entered, index is always 0. The randoms 172 // contains the last 12 numbers in the seque 165 // contains the last 12 numbers in the sequents: s[0] is x[a+11], 173 // s[1] is x[a+10] ... and s[11] is x[a] for 166 // s[1] is x[a+10] ... and s[11] is x[a] for some a. Carry contains 174 // the last carry value (c[a+11]). 167 // the last carry value (c[a+11]). 175 // 168 // 176 // The recursion relation (3) in Luscher's n 169 // The recursion relation (3) in Luscher's note says 177 // delta[n] = x[n-s] = x[n-r] -c[n-1] or f 170 // delta[n] = x[n-s] = x[n-r] -c[n-1] or for n=a+12, 178 // delta[a+12] = x[a+7] - x[a] -c[a+11] wh 171 // delta[a+12] = x[a+7] - x[a] -c[a+11] where we use r=12, s=5 per eqn. (7) 179 // This reduces to 172 // This reduces to 180 // s[11] = s[4] - s[11] - carry. 173 // s[11] = s[4] - s[11] - carry. 181 // The next number similarly will be given b 174 // The next number similarly will be given by s[10] = s[3] - s[10] - carry, 182 // and so forth until s[0] is filled. 175 // and so forth until s[0] is filled. 183 // 176 // 184 // However, we need to skip 397, 202 or 109 177 // However, we need to skip 397, 202 or 109 numbers - these are not divisible 185 // by 12 - to "fare well in the spectral tes 178 // by 12 - to "fare well in the spectral test". 186 179 187 advance(pDozens); 180 advance(pDozens); 188 181 189 // Since we wish at the end to have the 12 l 182 // Since we wish at the end to have the 12 last numbers in the order of 190 // s[11] first, till s[0] last, we will have 183 // s[11] first, till s[0] last, we will have to do 1, 10, or 1 iterations 191 // and then re-arrange to place to get the o 184 // and then re-arrange to place to get the oldest one in s[11]. 192 // Generically, this will imply re-arranging 185 // Generically, this will imply re-arranging the s array at the end, 193 // but we can treat the special case of endI 186 // but we can treat the special case of endIters = 1 separately for superior 194 // efficiency in the cases of levels 0 and 2 187 // efficiency in the cases of levels 0 and 2. 195 188 196 double y1; << 189 register double y1; 197 190 198 if ( endIters == 1 ) { // Luxury levels 0 191 if ( endIters == 1 ) { // Luxury levels 0 and 2 will go here 199 y1 = randoms[ 4] - randoms[11] - carry; 192 y1 = randoms[ 4] - randoms[11] - carry; 200 if ( y1 < 0.0 ) { 193 if ( y1 < 0.0 ) { 201 y1 += 1.0; 194 y1 += 1.0; 202 carry = twoToMinus_48(); 195 carry = twoToMinus_48(); 203 } else { 196 } else { 204 carry = 0.0; 197 carry = 0.0; 205 } 198 } 206 randoms[11] = randoms[10]; 199 randoms[11] = randoms[10]; 207 randoms[10] = randoms[ 9]; 200 randoms[10] = randoms[ 9]; 208 randoms[ 9] = randoms[ 8]; 201 randoms[ 9] = randoms[ 8]; 209 randoms[ 8] = randoms[ 7]; 202 randoms[ 8] = randoms[ 7]; 210 randoms[ 7] = randoms[ 6]; 203 randoms[ 7] = randoms[ 6]; 211 randoms[ 6] = randoms[ 5]; 204 randoms[ 6] = randoms[ 5]; 212 randoms[ 5] = randoms[ 4]; 205 randoms[ 5] = randoms[ 4]; 213 randoms[ 4] = randoms[ 3]; 206 randoms[ 4] = randoms[ 3]; 214 randoms[ 3] = randoms[ 2]; 207 randoms[ 3] = randoms[ 2]; 215 randoms[ 2] = randoms[ 1]; 208 randoms[ 2] = randoms[ 1]; 216 randoms[ 1] = randoms[ 0]; 209 randoms[ 1] = randoms[ 0]; 217 randoms[ 0] = y1; 210 randoms[ 0] = y1; 218 211 219 } else { 212 } else { 220 213 221 int m, nr, ns; 214 int m, nr, ns; 222 for ( m = 0, nr = 11, ns = 4; m < endIters 215 for ( m = 0, nr = 11, ns = 4; m < endIters; ++m, --nr ) { 223 y1 = randoms [ns] - randoms[nr] - carry; 216 y1 = randoms [ns] - randoms[nr] - carry; 224 if ( y1 < 0.0 ) { 217 if ( y1 < 0.0 ) { 225 y1 += 1.0; 218 y1 += 1.0; 226 carry = twoToMinus_48(); 219 carry = twoToMinus_48(); 227 } else { 220 } else { 228 carry = 0.0; 221 carry = 0.0; 229 } 222 } 230 randoms[nr] = y1; 223 randoms[nr] = y1; 231 --ns; 224 --ns; 232 if ( ns < 0 ) { 225 if ( ns < 0 ) { 233 ns = 11; 226 ns = 11; 234 } 227 } 235 } // loop on m 228 } // loop on m 236 229 237 double temp[12]; 230 double temp[12]; 238 for (m=0; m<12; m++) { 231 for (m=0; m<12; m++) { 239 temp[m]=randoms[m]; 232 temp[m]=randoms[m]; 240 } 233 } 241 234 242 ns = 11 - endIters; 235 ns = 11 - endIters; 243 for (m=11; m>=0; --m) { 236 for (m=11; m>=0; --m) { 244 randoms[m] = temp[ns]; 237 randoms[m] = temp[ns]; 245 --ns; 238 --ns; 246 if ( ns < 0 ) { 239 if ( ns < 0 ) { 247 ns = 11; 240 ns = 11; 248 } 241 } 249 } 242 } 250 243 251 } 244 } 252 245 253 // Now when we return, there are 12 fresh us 246 // Now when we return, there are 12 fresh usable numbers in s[11] ... s[0] 254 247 255 index = 12; << 248 index = 11; 256 249 257 } // update() 250 } // update() 258 251 259 void Ranlux64Engine::advance(int dozens) { 252 void Ranlux64Engine::advance(int dozens) { 260 253 261 double y1, y2, y3; << 254 register double y1, y2, y3; 262 double cValue = twoToMinus_48(); << 255 register double cValue = twoToMinus_48(); 263 double zero = 0.0; << 256 register double zero = 0.0; 264 double one = 1.0; << 257 register double one = 1.0; 265 258 266 // Technical note: We use Luscher's trick 259 // Technical note: We use Luscher's trick to only do the 267 // carry subtraction when we really have t 260 // carry subtraction when we really have to. Like him, we use 268 // three registers instead of two so that 261 // three registers instead of two so that we avoid sequences 269 // like storing y1 then immediately replac 262 // like storing y1 then immediately replacing its value: 270 // some architectures lose time when this 263 // some architectures lose time when this is done. 271 264 272 // Luscher's ranlxd.c fills the stash go 265 // Luscher's ranlxd.c fills the stash going 273 // upward. We fill it downward to save a 266 // upward. We fill it downward to save a bit of time in the 274 // flat() routine at no cost later. This 267 // flat() routine at no cost later. This means that while 275 // Luscher's ir is jr+5, our n-r is (n-s)- 268 // Luscher's ir is jr+5, our n-r is (n-s)-5. (Note that 276 // though ranlxd.c initializes ir and jr t 269 // though ranlxd.c initializes ir and jr to 11 and 7, ir as 277 // used is 5 more than jr because update i 270 // used is 5 more than jr because update is entered after 278 // incrementing ir.) 271 // incrementing ir.) 279 // 272 // 280 273 281 // I have CAREFULLY checked that the algor 274 // I have CAREFULLY checked that the algorithms do match 282 // in all details. 275 // in all details. 283 276 284 int k; 277 int k; 285 for ( k = dozens; k > 0; --k ) { 278 for ( k = dozens; k > 0; --k ) { 286 279 287 y1 = randoms[ 4] - randoms[11] - carry; 280 y1 = randoms[ 4] - randoms[11] - carry; 288 281 289 y2 = randoms[ 3] - randoms[10]; 282 y2 = randoms[ 3] - randoms[10]; 290 if ( y1 < zero ) { 283 if ( y1 < zero ) { 291 y1 += one; 284 y1 += one; 292 y2 -= cValue; 285 y2 -= cValue; 293 } 286 } 294 randoms[11] = y1; 287 randoms[11] = y1; 295 288 296 y3 = randoms[ 2] - randoms[ 9]; 289 y3 = randoms[ 2] - randoms[ 9]; 297 if ( y2 < zero ) { 290 if ( y2 < zero ) { 298 y2 += one; 291 y2 += one; 299 y3 -= cValue; 292 y3 -= cValue; 300 } 293 } 301 randoms[10] = y2; 294 randoms[10] = y2; 302 295 303 y1 = randoms[ 1] - randoms[ 8]; 296 y1 = randoms[ 1] - randoms[ 8]; 304 if ( y3 < zero ) { 297 if ( y3 < zero ) { 305 y3 += one; 298 y3 += one; 306 y1 -= cValue; 299 y1 -= cValue; 307 } 300 } 308 randoms[ 9] = y3; 301 randoms[ 9] = y3; 309 302 310 y2 = randoms[ 0] - randoms[ 7]; 303 y2 = randoms[ 0] - randoms[ 7]; 311 if ( y1 < zero ) { 304 if ( y1 < zero ) { 312 y1 += one; 305 y1 += one; 313 y2 -= cValue; 306 y2 -= cValue; 314 } 307 } 315 randoms[ 8] = y1; 308 randoms[ 8] = y1; 316 309 317 y3 = randoms[11] - randoms[ 6]; 310 y3 = randoms[11] - randoms[ 6]; 318 if ( y2 < zero ) { 311 if ( y2 < zero ) { 319 y2 += one; 312 y2 += one; 320 y3 -= cValue; 313 y3 -= cValue; 321 } 314 } 322 randoms[ 7] = y2; 315 randoms[ 7] = y2; 323 316 324 y1 = randoms[10] - randoms[ 5]; 317 y1 = randoms[10] - randoms[ 5]; 325 if ( y3 < zero ) { 318 if ( y3 < zero ) { 326 y3 += one; 319 y3 += one; 327 y1 -= cValue; 320 y1 -= cValue; 328 } 321 } 329 randoms[ 6] = y3; 322 randoms[ 6] = y3; 330 323 331 y2 = randoms[ 9] - randoms[ 4]; 324 y2 = randoms[ 9] - randoms[ 4]; 332 if ( y1 < zero ) { 325 if ( y1 < zero ) { 333 y1 += one; 326 y1 += one; 334 y2 -= cValue; 327 y2 -= cValue; 335 } 328 } 336 randoms[ 5] = y1; 329 randoms[ 5] = y1; 337 330 338 y3 = randoms[ 8] - randoms[ 3]; 331 y3 = randoms[ 8] - randoms[ 3]; 339 if ( y2 < zero ) { 332 if ( y2 < zero ) { 340 y2 += one; 333 y2 += one; 341 y3 -= cValue; 334 y3 -= cValue; 342 } 335 } 343 randoms[ 4] = y2; 336 randoms[ 4] = y2; 344 337 345 y1 = randoms[ 7] - randoms[ 2]; 338 y1 = randoms[ 7] - randoms[ 2]; 346 if ( y3 < zero ) { 339 if ( y3 < zero ) { 347 y3 += one; 340 y3 += one; 348 y1 -= cValue; 341 y1 -= cValue; 349 } 342 } 350 randoms[ 3] = y3; 343 randoms[ 3] = y3; 351 344 352 y2 = randoms[ 6] - randoms[ 1]; 345 y2 = randoms[ 6] - randoms[ 1]; 353 if ( y1 < zero ) { 346 if ( y1 < zero ) { 354 y1 += one; 347 y1 += one; 355 y2 -= cValue; 348 y2 -= cValue; 356 } 349 } 357 randoms[ 2] = y1; 350 randoms[ 2] = y1; 358 351 359 y3 = randoms[ 5] - randoms[ 0]; 352 y3 = randoms[ 5] - randoms[ 0]; 360 if ( y2 < zero ) { 353 if ( y2 < zero ) { 361 y2 += one; 354 y2 += one; 362 y3 -= cValue; 355 y3 -= cValue; 363 } 356 } 364 randoms[ 1] = y2; 357 randoms[ 1] = y2; 365 358 366 if ( y3 < zero ) { 359 if ( y3 < zero ) { 367 y3 += one; 360 y3 += one; 368 carry = cValue; 361 carry = cValue; 369 } 362 } 370 randoms[ 0] = y3; 363 randoms[ 0] = y3; 371 364 372 } // End of major k loop doing 12 numbers at 365 } // End of major k loop doing 12 numbers at each cycle 373 366 374 } // advance(dozens) 367 } // advance(dozens) 375 368 376 void Ranlux64Engine::flatArray(const int size, 369 void Ranlux64Engine::flatArray(const int size, double* vect) { 377 for( int i=0; i < size; ++i ) { 370 for( int i=0; i < size; ++i ) { 378 vect[i] = flat(); 371 vect[i] = flat(); 379 } 372 } 380 } 373 } 381 374 382 void Ranlux64Engine::setSeed(long seed, int lu 375 void Ranlux64Engine::setSeed(long seed, int lux) { 383 376 384 // The initialization is carried out using a M 377 // The initialization is carried out using a Multiplicative 385 // Congruential generator using formula consta 378 // Congruential generator using formula constants of L'Ecuyer 386 // as described in "A review of pseudorandom n 379 // as described in "A review of pseudorandom number generators" 387 // (Fred James) published in Computer Physics 380 // (Fred James) published in Computer Physics Communications 60 (1990) 388 // pages 329-344 381 // pages 329-344 389 382 390 const int ecuyer_a(53668); 383 const int ecuyer_a(53668); 391 const int ecuyer_b(40014); 384 const int ecuyer_b(40014); 392 const int ecuyer_c(12211); 385 const int ecuyer_c(12211); 393 const int ecuyer_d(2147483563); 386 const int ecuyer_d(2147483563); 394 387 395 const int lux_levels[3] = {109, 202, 397}; 388 const int lux_levels[3] = {109, 202, 397}; 396 theSeed = seed; 389 theSeed = seed; 397 390 398 if( (lux > 2)||(lux < 0) ){ 391 if( (lux > 2)||(lux < 0) ){ 399 pDiscard = (lux >= 12) ? (lux-12) : lux_l 392 pDiscard = (lux >= 12) ? (lux-12) : lux_levels[1]; 400 }else{ 393 }else{ 401 pDiscard = lux_levels[luxury]; 394 pDiscard = lux_levels[luxury]; 402 } 395 } 403 pDozens = pDiscard / 12; 396 pDozens = pDiscard / 12; 404 endIters = pDiscard % 12; 397 endIters = pDiscard % 12; 405 398 406 long init_table[24]; 399 long init_table[24]; 407 long next_seed = seed; 400 long next_seed = seed; 408 long k_multiple; 401 long k_multiple; 409 int i; 402 int i; 410 next_seed &= 0xffffffff; 403 next_seed &= 0xffffffff; 411 while( next_seed >= ecuyer_d ) { 404 while( next_seed >= ecuyer_d ) { 412 next_seed -= ecuyer_d; 405 next_seed -= ecuyer_d; 413 } 406 } 414 407 415 for(i = 0;i != 24;i++){ 408 for(i = 0;i != 24;i++){ 416 k_multiple = next_seed / ecuyer_a; 409 k_multiple = next_seed / ecuyer_a; 417 next_seed = ecuyer_b * (next_seed - k_mul 410 next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a) 418 - k_mul 411 - k_multiple * ecuyer_c; 419 if(next_seed < 0) { 412 if(next_seed < 0) { 420 next_seed += ecuyer_d; 413 next_seed += ecuyer_d; 421 } 414 } 422 next_seed &= 0xffffffff; 415 next_seed &= 0xffffffff; 423 init_table[i] = next_seed; 416 init_table[i] = next_seed; 424 } 417 } 425 // are we on a 64bit machine? 418 // are we on a 64bit machine? 426 if( sizeof(long) >= 8 ) { 419 if( sizeof(long) >= 8 ) { 427 int64_t topbits1, topbits2; << 420 long topbits1, topbits2; 428 #ifdef WIN32 421 #ifdef WIN32 429 topbits1 = ( (int64_t) seed >> 32) & 0xff << 422 topbits1 = ( seed >> 32) & 0xffff ; 430 topbits2 = ( (int64_t) seed >> 48) & 0xff << 423 topbits2 = ( seed >> 48) & 0xffff ; 431 #else 424 #else 432 topbits1 = detail::rshift<32>(seed) & 0xf 425 topbits1 = detail::rshift<32>(seed) & 0xffff ; 433 topbits2 = detail::rshift<48>(seed) & 0xf 426 topbits2 = detail::rshift<48>(seed) & 0xffff ; 434 #endif 427 #endif 435 init_table[0] ^= topbits1; 428 init_table[0] ^= topbits1; 436 init_table[2] ^= topbits2; 429 init_table[2] ^= topbits2; 437 //std::cout << " init_table[0] " << init_ 430 //std::cout << " init_table[0] " << init_table[0] << " from " << topbits1 << std::endl; 438 //std::cout << " init_table[2] " << init_ 431 //std::cout << " init_table[2] " << init_table[2] << " from " << topbits2 << std::endl; 439 } 432 } 440 433 441 for(i = 0;i < 12; i++){ 434 for(i = 0;i < 12; i++){ 442 randoms[i] = (init_table[2*i ] ) * 435 randoms[i] = (init_table[2*i ] ) * 2.0 * twoToMinus_32() + 443 (init_table[2*i+1] >> 15) * 436 (init_table[2*i+1] >> 15) * twoToMinus_48(); 444 //if( randoms[i] < 0. || randoms[i] > 1. 437 //if( randoms[i] < 0. || randoms[i] > 1. ) { 445 //std::cout << "setSeed: init_table " << 438 //std::cout << "setSeed: init_table " << init_table[2*i ] << std::endl; 446 //std::cout << "setSeed: init_table " << 439 //std::cout << "setSeed: init_table " << init_table[2*i+1] << std::endl; 447 //std::cout << "setSeed: random " << i < 440 //std::cout << "setSeed: random " << i << " is " << randoms[i] << std::endl; 448 //} 441 //} 449 } 442 } 450 443 451 carry = 0.0; 444 carry = 0.0; 452 if ( randoms[11] == 0. ) carry = twoToMinus_ 445 if ( randoms[11] == 0. ) carry = twoToMinus_48(); 453 // Perform an update before returning the fi << 446 index = 11; 454 index = -1; << 455 447 456 } // setSeed() 448 } // setSeed() 457 449 458 void Ranlux64Engine::setSeeds(const long * see 450 void Ranlux64Engine::setSeeds(const long * seeds, int lux) { 459 // old code only uses the first long in seeds 451 // old code only uses the first long in seeds 460 // setSeed( *seeds ? *seeds : 32767, lux ); 452 // setSeed( *seeds ? *seeds : 32767, lux ); 461 // theSeeds = seeds; 453 // theSeeds = seeds; 462 454 463 // using code from Ranlux - even those are 32b 455 // using code from Ranlux - even those are 32bit seeds, 464 // that is good enough to completely different 456 // that is good enough to completely differentiate the sequences 465 457 466 const int ecuyer_a = 53668; 458 const int ecuyer_a = 53668; 467 const int ecuyer_b = 40014; 459 const int ecuyer_b = 40014; 468 const int ecuyer_c = 12211; 460 const int ecuyer_c = 12211; 469 const int ecuyer_d = 2147483563; 461 const int ecuyer_d = 2147483563; 470 462 471 const int lux_levels[3] = {109, 202, 397}; 463 const int lux_levels[3] = {109, 202, 397}; 472 const long *seedptr; 464 const long *seedptr; 473 465 474 theSeeds = seeds; 466 theSeeds = seeds; 475 seedptr = seeds; 467 seedptr = seeds; 476 468 477 if(seeds == 0){ 469 if(seeds == 0){ 478 setSeed(theSeed,lux); 470 setSeed(theSeed,lux); 479 theSeeds = &theSeed; 471 theSeeds = &theSeed; 480 return; 472 return; 481 } 473 } 482 474 483 theSeed = *seeds; 475 theSeed = *seeds; 484 476 485 // number of additional random numbers that ne 477 // number of additional random numbers that need to be 'thrown away' 486 // every 24 numbers is set using luxury level 478 // every 24 numbers is set using luxury level variable. 487 479 488 if( (lux > 2)||(lux < 0) ){ 480 if( (lux > 2)||(lux < 0) ){ 489 pDiscard = (lux >= 12) ? (lux-12) : lux_l 481 pDiscard = (lux >= 12) ? (lux-12) : lux_levels[1]; 490 }else{ 482 }else{ 491 pDiscard = lux_levels[luxury]; 483 pDiscard = lux_levels[luxury]; 492 } 484 } 493 pDozens = pDiscard / 12; 485 pDozens = pDiscard / 12; 494 endIters = pDiscard % 12; 486 endIters = pDiscard % 12; 495 487 496 long init_table[24]; 488 long init_table[24]; 497 long next_seed = *seeds; 489 long next_seed = *seeds; 498 long k_multiple; 490 long k_multiple; 499 int i; 491 int i; 500 492 501 for( i = 0;(i != 24)&&(*seedptr != 0);i++){ 493 for( i = 0;(i != 24)&&(*seedptr != 0);i++){ 502 init_table[i] = *seedptr & 0xffffffff; 494 init_table[i] = *seedptr & 0xffffffff; 503 seedptr++; 495 seedptr++; 504 } 496 } 505 497 506 if(i != 24){ 498 if(i != 24){ 507 next_seed = init_table[i-1]; 499 next_seed = init_table[i-1]; 508 for(;i != 24;i++){ 500 for(;i != 24;i++){ 509 k_multiple = next_seed / ecuyer_a; 501 k_multiple = next_seed / ecuyer_a; 510 next_seed = ecuyer_b * (next_seed - k_multip 502 next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a) 511 - k_multip 503 - k_multiple * ecuyer_c; 512 if(next_seed < 0) { 504 if(next_seed < 0) { 513 next_seed += ecuyer_d; 505 next_seed += ecuyer_d; 514 } 506 } 515 next_seed &= 0xffffffff; 507 next_seed &= 0xffffffff; 516 init_table[i] = next_seed; 508 init_table[i] = next_seed; 517 } 509 } 518 } 510 } 519 511 520 for(i = 0;i < 12; i++){ 512 for(i = 0;i < 12; i++){ 521 randoms[i] = (init_table[2*i ] ) * 513 randoms[i] = (init_table[2*i ] ) * 2.0 * twoToMinus_32() + 522 (init_table[2*i+1] >> 15) * 514 (init_table[2*i+1] >> 15) * twoToMinus_48(); 523 } 515 } 524 516 525 carry = 0.0; 517 carry = 0.0; 526 if ( randoms[11] == 0. ) carry = twoToMinus_ 518 if ( randoms[11] == 0. ) carry = twoToMinus_48(); 527 // Perform an update before returning the fi << 519 index = 11; 528 index = -1; << 529 520 530 } 521 } 531 522 532 void Ranlux64Engine::saveStatus( const char fi 523 void Ranlux64Engine::saveStatus( const char filename[] ) const 533 { 524 { 534 std::ofstream outFile( filename, std::ios:: 525 std::ofstream outFile( filename, std::ios::out ) ; 535 if (!outFile.bad()) { 526 if (!outFile.bad()) { 536 outFile << "Uvec\n"; 527 outFile << "Uvec\n"; 537 std::vector<unsigned long> v = put(); 528 std::vector<unsigned long> v = put(); 538 for (unsigned int i=0; i<v.size(); ++i) { 529 for (unsigned int i=0; i<v.size(); ++i) { 539 outFile << v[i] << "\n"; 530 outFile << v[i] << "\n"; 540 } 531 } 541 } 532 } 542 } 533 } 543 534 544 void Ranlux64Engine::restoreStatus( const char 535 void Ranlux64Engine::restoreStatus( const char filename[] ) 545 { 536 { 546 std::ifstream inFile( filename, std::ios::i 537 std::ifstream inFile( filename, std::ios::in); 547 if (!checkFile ( inFile, filename, engineNa 538 if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) { 548 std::cerr << " -- Engine state remains u 539 std::cerr << " -- Engine state remains unchanged\n"; 549 return; 540 return; 550 } 541 } 551 if ( possibleKeywordInput ( inFile, "Uvec", 542 if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) { 552 std::vector<unsigned long> v; 543 std::vector<unsigned long> v; 553 unsigned long xin; 544 unsigned long xin; 554 for (unsigned int ivec=0; ivec < VECTOR_ST 545 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) { 555 inFile >> xin; 546 inFile >> xin; 556 if (!inFile) { 547 if (!inFile) { 557 inFile.clear(std::ios::badbit | inFile 548 inFile.clear(std::ios::badbit | inFile.rdstate()); 558 std::cerr << "\nJamesRandom state (vec 549 std::cerr << "\nJamesRandom state (vector) description improper." 559 << "\nrestoreStatus has failed." 550 << "\nrestoreStatus has failed." 560 << "\nInput stream is probably mispos 551 << "\nInput stream is probably mispositioned now." << std::endl; 561 return; 552 return; 562 } 553 } 563 v.push_back(xin); 554 v.push_back(xin); 564 } 555 } 565 getState(v); 556 getState(v); 566 return; 557 return; 567 } 558 } 568 559 569 if (!inFile.bad() && !inFile.eof()) { 560 if (!inFile.bad() && !inFile.eof()) { 570 // inFile >> theSeed; removed -- encompas 561 // inFile >> theSeed; removed -- encompased by possibleKeywordInput 571 for (int i=0; i<12; ++i) { 562 for (int i=0; i<12; ++i) { 572 inFile >> randoms[i]; 563 inFile >> randoms[i]; 573 } 564 } 574 inFile >> carry; inFile >> index; 565 inFile >> carry; inFile >> index; 575 inFile >> luxury; inFile >> pDiscard; 566 inFile >> luxury; inFile >> pDiscard; 576 pDozens = pDiscard / 12; 567 pDozens = pDiscard / 12; 577 endIters = pDiscard % 12; 568 endIters = pDiscard % 12; 578 } 569 } 579 } 570 } 580 571 581 void Ranlux64Engine::showStatus() const 572 void Ranlux64Engine::showStatus() const 582 { 573 { 583 std::cout << std::endl; 574 std::cout << std::endl; 584 std::cout << "--------- Ranlux engine statu 575 std::cout << "--------- Ranlux engine status ---------" << std::endl; 585 std::cout << " Initial seed = " << theSeed 576 std::cout << " Initial seed = " << theSeed << std::endl; 586 std::cout << " randoms[] = "; 577 std::cout << " randoms[] = "; 587 for (int i=0; i<12; ++i) { 578 for (int i=0; i<12; ++i) { 588 std::cout << randoms[i] << std::endl; 579 std::cout << randoms[i] << std::endl; 589 } 580 } 590 std::cout << std::endl; 581 std::cout << std::endl; 591 std::cout << " carry = " << carry << ", ind 582 std::cout << " carry = " << carry << ", index = " << index << std::endl; 592 std::cout << " luxury = " << luxury << " pD 583 std::cout << " luxury = " << luxury << " pDiscard = " 593 << pDiscard << std::endl; 584 << pDiscard << std::endl; 594 std::cout << "----------------------------- 585 std::cout << "----------------------------------------" << std::endl; 595 } 586 } 596 587 597 std::ostream & Ranlux64Engine::put( std::ostre 588 std::ostream & Ranlux64Engine::put( std::ostream& os ) const 598 { 589 { 599 char beginMarker[] = "Ranlux64Engine-begin" 590 char beginMarker[] = "Ranlux64Engine-begin"; 600 os << beginMarker << "\nUvec\n"; 591 os << beginMarker << "\nUvec\n"; 601 std::vector<unsigned long> v = put(); 592 std::vector<unsigned long> v = put(); 602 for (unsigned int i=0; i<v.size(); ++i) { 593 for (unsigned int i=0; i<v.size(); ++i) { 603 os << v[i] << "\n"; 594 os << v[i] << "\n"; 604 } 595 } 605 return os; 596 return os; 606 } 597 } 607 598 608 std::vector<unsigned long> Ranlux64Engine::put 599 std::vector<unsigned long> Ranlux64Engine::put () const { 609 std::vector<unsigned long> v; 600 std::vector<unsigned long> v; 610 v.push_back (engineIDulong<Ranlux64Engine>() 601 v.push_back (engineIDulong<Ranlux64Engine>()); 611 std::vector<unsigned long> t; 602 std::vector<unsigned long> t; 612 for (int i=0; i<12; ++i) { 603 for (int i=0; i<12; ++i) { 613 t = DoubConv::dto2longs(randoms[i]); 604 t = DoubConv::dto2longs(randoms[i]); 614 v.push_back(t[0]); v.push_back(t[1]); 605 v.push_back(t[0]); v.push_back(t[1]); 615 } 606 } 616 t = DoubConv::dto2longs(carry); 607 t = DoubConv::dto2longs(carry); 617 v.push_back(t[0]); v.push_back(t[1]); 608 v.push_back(t[0]); v.push_back(t[1]); 618 v.push_back(static_cast<unsigned long>(index 609 v.push_back(static_cast<unsigned long>(index)); 619 v.push_back(static_cast<unsigned long>(luxur 610 v.push_back(static_cast<unsigned long>(luxury)); 620 v.push_back(static_cast<unsigned long>(pDisc 611 v.push_back(static_cast<unsigned long>(pDiscard)); 621 return v; 612 return v; 622 } 613 } 623 614 624 std::istream & Ranlux64Engine::get ( std::istr 615 std::istream & Ranlux64Engine::get ( std::istream& is ) 625 { 616 { 626 char beginMarker [MarkerLen]; 617 char beginMarker [MarkerLen]; 627 is >> std::ws; 618 is >> std::ws; 628 is.width(MarkerLen); // causes the next rea 619 is.width(MarkerLen); // causes the next read to the char* to be <= 629 // that many bytes, INCLUDING A TERMINAT 620 // that many bytes, INCLUDING A TERMINATION \0 630 // (Stroustrup, section 21.3.2) 621 // (Stroustrup, section 21.3.2) 631 is >> beginMarker; 622 is >> beginMarker; 632 if (strcmp(beginMarker,"Ranlux64Engine-begin 623 if (strcmp(beginMarker,"Ranlux64Engine-begin")) { 633 is.clear(std::ios::badbit | is.rdstate()) 624 is.clear(std::ios::badbit | is.rdstate()); 634 std::cerr << "\nInput stream mispositione 625 std::cerr << "\nInput stream mispositioned or" 635 << "\nRanlux64Engine state descriptio 626 << "\nRanlux64Engine state description missing or" 636 << "\nwrong engine type found." << st 627 << "\nwrong engine type found." << std::endl; 637 return is; 628 return is; 638 } 629 } 639 return getState(is); 630 return getState(is); 640 } 631 } 641 632 642 std::string Ranlux64Engine::beginTag ( ) { 633 std::string Ranlux64Engine::beginTag ( ) { 643 return "Ranlux64Engine-begin"; 634 return "Ranlux64Engine-begin"; 644 } 635 } 645 636 646 std::istream & Ranlux64Engine::getState ( std: 637 std::istream & Ranlux64Engine::getState ( std::istream& is ) 647 { 638 { 648 if ( possibleKeywordInput ( is, "Uvec", theS 639 if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) { 649 std::vector<unsigned long> v; 640 std::vector<unsigned long> v; 650 unsigned long uu; 641 unsigned long uu; 651 for (unsigned int ivec=0; ivec < VECTOR_ST 642 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) { 652 is >> uu; 643 is >> uu; 653 if (!is) { 644 if (!is) { 654 is.clear(std::ios::badbit | is.rdstate 645 is.clear(std::ios::badbit | is.rdstate()); 655 std::cerr << "\nRanlux64Engine state ( 646 std::cerr << "\nRanlux64Engine state (vector) description improper." 656 << "\ngetState() has failed." 647 << "\ngetState() has failed." 657 << "\nInput stream is probably mispos 648 << "\nInput stream is probably mispositioned now." << std::endl; 658 return is; 649 return is; 659 } 650 } 660 v.push_back(uu); 651 v.push_back(uu); 661 } 652 } 662 getState(v); 653 getState(v); 663 return (is); 654 return (is); 664 } 655 } 665 656 666 // is >> theSeed; Removed, encompassed by po 657 // is >> theSeed; Removed, encompassed by possibleKeywordInput() 667 658 668 char endMarker [MarkerLen]; 659 char endMarker [MarkerLen]; 669 for (int i=0; i<12; ++i) { 660 for (int i=0; i<12; ++i) { 670 is >> randoms[i]; 661 is >> randoms[i]; 671 } 662 } 672 is >> carry; is >> index; 663 is >> carry; is >> index; 673 is >> luxury; is >> pDiscard; 664 is >> luxury; is >> pDiscard; 674 pDozens = pDiscard / 12; 665 pDozens = pDiscard / 12; 675 endIters = pDiscard % 12; 666 endIters = pDiscard % 12; 676 is >> std::ws; 667 is >> std::ws; 677 is.width(MarkerLen); 668 is.width(MarkerLen); 678 is >> endMarker; 669 is >> endMarker; 679 if (strcmp(endMarker,"Ranlux64Engine-end")) 670 if (strcmp(endMarker,"Ranlux64Engine-end")) { 680 is.clear(std::ios::badbit | is.rdstate()) 671 is.clear(std::ios::badbit | is.rdstate()); 681 std::cerr << "\nRanlux64Engine state desc 672 std::cerr << "\nRanlux64Engine state description incomplete." 682 << "\nInput stream is probably mispos 673 << "\nInput stream is probably mispositioned now." << std::endl; 683 return is; 674 return is; 684 } 675 } 685 return is; 676 return is; 686 } 677 } 687 678 688 bool Ranlux64Engine::get (const std::vector<un 679 bool Ranlux64Engine::get (const std::vector<unsigned long> & v) { 689 if ((v[0] & 0xffffffffUL) != engineIDulong<R 680 if ((v[0] & 0xffffffffUL) != engineIDulong<Ranlux64Engine>()) { 690 std::cerr << 681 std::cerr << 691 "\nRanlux64Engine get:state vector has w 682 "\nRanlux64Engine get:state vector has wrong ID word - state unchanged\n"; 692 return false; 683 return false; 693 } 684 } 694 return getState(v); 685 return getState(v); 695 } 686 } 696 687 697 bool Ranlux64Engine::getState (const std::vect 688 bool Ranlux64Engine::getState (const std::vector<unsigned long> & v) { 698 if (v.size() != VECTOR_STATE_SIZE ) { 689 if (v.size() != VECTOR_STATE_SIZE ) { 699 std::cerr << 690 std::cerr << 700 "\nRanlux64Engine get:state vector has w 691 "\nRanlux64Engine get:state vector has wrong length - state unchanged\n"; 701 return false; 692 return false; 702 } 693 } 703 std::vector<unsigned long> t(2); 694 std::vector<unsigned long> t(2); 704 for (int i=0; i<12; ++i) { 695 for (int i=0; i<12; ++i) { 705 t[0] = v[2*i+1]; t[1] = v[2*i+2]; 696 t[0] = v[2*i+1]; t[1] = v[2*i+2]; 706 randoms[i] = DoubConv::longs2double(t); 697 randoms[i] = DoubConv::longs2double(t); 707 } 698 } 708 t[0] = v[25]; t[1] = v[26]; 699 t[0] = v[25]; t[1] = v[26]; 709 carry = DoubConv::longs2double(t); 700 carry = DoubConv::longs2double(t); 710 index = (int)v[27]; << 701 index = v[27]; 711 luxury = (int)v[28]; << 702 luxury = v[28]; 712 pDiscard = (int)v[29]; << 703 pDiscard = v[29]; 713 return true; 704 return true; 714 } 705 } 715 706 716 } // namespace CLHEP 707 } // namespace CLHEP 717 708