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