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