| Geant4 Cross Reference |
1 //
1
2 // -*- C++ -*-
3 //
4 // -------------------------------------------
5 // HEP Random
6 // --- MixMaxRng ---
7 // class implementation fi
8 // -------------------------------------------
9 //
10 // This file interfaces the MixMax PseudoRando
11 // proposed by:
12 //
13 // G.K.Savvidy and N.G.Ter-Arutyunian,
14 // On the Monte Carlo simulation of physical
15 // J.Comput.Phys. 97, 566 (1991);
16 // Preprint EPI-865-16-86, Yerevan, Jan. 198
17 // http://dx.doi.org/10.1016/0021-9991(91)90
18 //
19 // K.Savvidy
20 // "The MIXMAX random number generator"
21 // Comp. Phys. Commun. (2015)
22 // http://dx.doi.org/10.1016/j.cpc.2015.06.0
23 //
24 // K.Savvidy and G.Savvidy
25 // "Spectrum and Entropy of C-systems. MIXMA
26 // Chaos, Solitons & Fractals, Volume 91, (2
27 // http://dx.doi.org/10.1016/j.chaos.2016.05
28 //
29 // ===========================================
30 // Implementation by Konstantin Savvidy - Copy
31 // July 2023 - Updated class structure upon su
32 // September 2023 - fix (re-)initialization fr
33 // ===========================================
34
35 #include "CLHEP/Random/Random.h"
36 #include "CLHEP/Random/MixMaxRng.h"
37 #include "CLHEP/Random/engineIDulong.h"
38 #include "CLHEP/Utility/atomic_int.h"
39
40 #include <atomic>
41 #include <cmath>
42 #include <algorithm>
43 #include <iostream>
44 #include <string.h> // for strcmp
45 #include <vector>
46
47 const unsigned long MASK32=0xffffffff;
48
49 namespace CLHEP {
50
51 namespace {
52 // Number of instances with automatic seed s
53 CLHEP_ATOMIC_INT_TYPE numberOfEngines(0);
54 }
55
56 static const int MarkerLen = 64; // Enough roo
57
58 MixMaxRng::MixMaxRng()
59 : HepRandomEngine()
60 {
61 int numEngines = ++numberOfEngines;
62 setSeed(static_cast<long>(numEngines));
63 }
64
65 MixMaxRng::MixMaxRng(long seed)
66 : HepRandomEngine()
67 {
68 theSeed=seed;
69 setSeed(seed);
70 }
71
72 MixMaxRng::MixMaxRng(std::istream& is)
73 : HepRandomEngine()
74 {
75 get(is);
76 }
77
78 MixMaxRng::~MixMaxRng()
79 {
80 }
81
82 MixMaxRng::MixMaxRng(const MixMaxRng& rng)
83 : HepRandomEngine(rng)
84 {
85 std::copy(rng.V, rng.V+N, V);
86 sumtot= rng.sumtot;
87 counter= rng.counter;
88 }
89
90 MixMaxRng& MixMaxRng::operator=(const MixMaxRn
91 {
92 // Check assignment to self
93 //
94 if (this == &rng) { return *this; }
95
96 // Copy base class data
97 //
98 HepRandomEngine::operator=(rng);
99
100 std::copy(rng.V, rng.V+N, V);
101 sumtot= rng.sumtot;
102 counter= rng.counter;
103
104 return *this;
105 }
106
107 void MixMaxRng::saveStatus( const char filenam
108 {
109 // Create a C file-handle or an object that
110 FILE *fh= fopen(filename, "w");
111 if( fh )
112 {
113 int j;
114 fprintf(fh, "mixmax state, file version 1
115 fprintf(fh, "N=%u; V[N]={", rng_get_N() )
116 for (j=0; (j< (rng_get_N()-1) ); ++j) {
117 fprintf(fh, "%llu, ", (unsigned long
118 }
119 fprintf(fh, "%llu", (unsigned long long)V
120 fprintf(fh, "}; " );
121 fprintf(fh, "counter=%u; ", counter );
122 fprintf(fh, "sumtot=%llu;\n", (unsigned l
123 fclose(fh);
124 }
125 }
126
127 void MixMaxRng::restoreStatus( const char file
128 {
129 // a function for reading the state from a
130 FILE* fin;
131 if( ( fin = fopen(filename, "r") ) )
132 {
133 char l=0;
134 while ( l != '{' ) { // 0x7B = "{"
135 l=fgetc(fin); // proceed until hitting
136 }
137 ungetc(' ', fin);
138 }
139 else
140 {
141 fprintf(stderr, "mixmax -> read_state: e
142 throw std::runtime_error("Error in readi
143 }
144
145 myuint_t vecVal;
146 //printf("mixmax -> read_state: starting to
147 if (!fscanf(fin, "%llu", (unsigned long lon
148 {
149 fprintf(stderr, "mixmax -> read_state: er
150 throw std::runtime_error("Error in readin
151 }
152
153 for (int i = 1; i < rng_get_N(); ++i)
154 {
155 if (!fscanf(fin, ", %llu", (unsigned long
156 {
157 fprintf(stderr, "mixmax -> read_state:
158 throw std::runtime_error("Error in read
159 }
160 if( vecVal <= MixMaxRng::M61 )
161 {
162 V[i] = vecVal;
163 }
164 else
165 {
166 fprintf(stderr, "mixmax -> read_state:
167 " ( must be less than %llu ) "
168 " obtained from reading file %s
169 , (unsigned long long)vecVal, (
170 }
171 }
172
173 int incounter;
174 if (!fscanf( fin, "}; counter=%i; ", &incou
175 {
176 fprintf(stderr, "mixmax -> read_state: er
177 throw std::runtime_error("Error in readin
178 }
179 if( incounter <= rng_get_N() )
180 {
181 counter = incounter;
182 }
183 else
184 {
185 fprintf(stderr, "mixmax -> read_state: In
186 " Must be 0 <= counter < %u\n" ,
187 print_state();
188 throw std::runtime_error("Error in readin
189 }
190 precalc();
191 myuint_t insumtot;
192 if (!fscanf( fin, "sumtot=%llu\n", (unsigne
193 {
194 fprintf(stderr, "mixmax -> read_state: er
195 throw std::runtime_error("Error in readin
196 }
197
198 if (sumtot != insumtot)
199 {
200 fprintf(stderr, "mixmax -> checksum error
201 throw std::runtime_error("Error in readin
202 }
203 fclose(fin);
204 }
205
206 void MixMaxRng::showStatus() const
207 {
208 std::cout << std::endl;
209 std::cout << "------- MixMaxRng engine stat
210
211 std::cout << " Current state vector is:" <<
212 print_state();
213 std::cout << "-----------------------------
214 }
215
216 // Preferred Seeding method
217 // the values of 'Seeds' must be valid 32-bit
218 // Higher order bits will be ignored!!
219
220 void MixMaxRng::setSeeds(const long* Seeds, in
221 {
222 myID_t seed0, seed1= 0, seed2= 0, seed3= 0;
223
224 if( seedNum < 1 ) { // Assuming at least 2
225 seed0= static_cast<myID_t>(Seeds[0]) &
226 seed1= static_cast<myID_t>(Seeds[1]) &
227 }
228 else
229 {
230 if( seedNum < 4 ) {
231 seed0= static_cast<myID_t>(Seeds[0]) &
232 if( seedNum > 1){ seed1= static_cast<my
233 if( seedNum > 2){ seed2= static_cast<my
234 }
235 if( seedNum >= 4 ) {
236 seed0= static_cast<myID_t>(Seeds[0]) &
237 seed1= static_cast<myID_t>(Seeds[1]) &
238 seed2= static_cast<myID_t>(Seeds[2]) &
239 seed3= static_cast<myID_t>(Seeds[3]) &
240 }
241 }
242 theSeed = Seeds[0];
243 theSeeds = Seeds;
244 seed_uniquestream(seed3, seed2, seed1, seed
245 }
246
247 std::string MixMaxRng::engineName()
248 {
249 return "MixMaxRng";
250 }
251
252 constexpr int MixMaxRng::rng_get_N()
253 {
254 return N;
255 }
256
257 void MixMaxRng::flatArray(const int size, doub
258 {
259 for (int i=0; i<size; ++i) { vect[i] = flat
260 }
261
262 std::ostream & MixMaxRng::put ( std::ostream&
263 {
264 char beginMarker[] = "MixMaxRng-begin";
265 char endMarker[] = "MixMaxRng-end";
266
267 long pr = os.precision(24);
268 os << beginMarker << " ";
269 os << theSeed << "\n";
270 for (int i=0; i<rng_get_N(); ++i) {
271 os << V[i] << "\n";
272 }
273 os << counter << "\n";
274 os << sumtot << "\n";
275 os << endMarker << "\n";
276 os.precision(pr);
277 return os;
278 }
279
280 std::vector<unsigned long> MixMaxRng::put () c
281 {
282 std::vector<unsigned long> vec;
283 vec.push_back (engineIDulong<MixMaxRng>());
284 for (int i=0; i<rng_get_N(); ++i)
285 {
286 vec.push_back(static_cast<unsigned long>(
287 // little-ended order on all platforms
288 vec.push_back(static_cast<unsigned long>(
289 // pack uint64 into a data structure wh
290 }
291 vec.push_back(static_cast<unsigned long>(co
292 vec.push_back(static_cast<unsigned long>(su
293 vec.push_back(static_cast<unsigned long>(su
294 return vec;
295 }
296
297 std::istream & MixMaxRng::get ( std::istream&
298 {
299 char beginMarker [MarkerLen];
300 is >> std::ws;
301 is.width(MarkerLen); // causes the next re
302 // that many bytes, I
303 // (Stroustrup, secti
304 is >> beginMarker;
305 if (strcmp(beginMarker,"MixMaxRng-begin"))
306 is.clear(std::ios::badbit | is.rdstate()
307 std::cerr << "\nInput stream misposition
308 << "\nMixMaxRng state descript
309 << "\nwrong engine type found.
310 return is;
311 }
312 return getState(is);
313 }
314
315 std::string MixMaxRng::beginTag ()
316 {
317 return "MixMaxRng-begin";
318 }
319
320 std::istream & MixMaxRng::getState ( std::ist
321 {
322 char endMarker[MarkerLen];
323 is >> theSeed;
324 for (int i=0; i<rng_get_N(); ++i) is >> V[
325 is >> counter;
326 myuint_t checksum;
327 is >> checksum;
328 is >> std::ws;
329 is.width(MarkerLen);
330 is >> endMarker;
331 if (strcmp(endMarker,"MixMaxRng-end")) {
332 is.clear(std::ios::badbit | is.rdstate(
333 std::cerr << "\nMixMaxRng state descrip
334 << "\nInput stream is probabl
335 return is;
336 }
337 if ( counter < 0 || counter > rng_get_N() )
338 std::cerr << "\nMixMaxRng::getState():
339 << "vector read wrong value o
340 << "\nInput stream is probabl
341 return is;
342 }
343 precalc();
344 if ( checksum != sumtot) {
345 std::cerr << "\nMixMaxRng::getState():
346 << "checksum disagrees with v
347 << "\nInput stream is probabl
348 return is;
349 }
350 return is;
351 }
352
353 bool MixMaxRng::get (const std::vector<unsigne
354 {
355 if ((vec[0] & 0xffffffffUL) != engineIDulon
356 {
357 std::cerr <<
358 "\nMixMaxRng::get(): vector has wrong
359 return false;
360 }
361 return getState(vec);
362 }
363
364 bool MixMaxRng::getState (const std::vector<un
365 {
366 if (vec.size() != VECTOR_STATE_SIZE ) {
367 std::cerr <<
368 "\nMixMaxRng::getState(): vector has w
369 return false;
370 }
371 for (int i=1; i<2*rng_get_N() ; i=i+2) {
372 V[i/2]= ( (vec[i] & MASK32) | ( (myuint_t
373 // unpack from a data structure which is
374 }
375 counter = (int)vec[2*rng_get_N()+1];
376 precalc();
377 if ( ( (vec[2*rng_get_N()+2] & MASK32)
378 | ( (myuint_t)(vec[2*rng_get_N()+3]) <
379 std::cerr << "\nMixMaxRng::getState(): ve
380 << "\nInput vector is probably
381 return false;
382 }
383 return true;
384 }
385
386 myuint_t MixMaxRng::iterate_raw_vec(myuint_t*
387 {
388 // operates with a raw vector, uses known s
389
390 myuint_t tempP, tempV;
391 Y[0] = ( tempV = sumtotOld);
392 myuint_t insumtot = Y[0], ovflow = 0; // wi
393 tempP = 0; // will keep a part
394 for (int i=1; (i<N); ++i)
395 {
396 myuint_t tempPO = MULWU(tempP);
397 tempP = modadd(tempP, Y[i]);
398 tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+t
399 Y[i] = tempV;
400 insumtot += tempV; if (insumtot < tempV)
401 }
402 return MIXMAX_MOD_MERSENNE(MIXMAX_MOD_MERSE
403 }
404
405 myuint_t MixMaxRng::get_next()
406 {
407 int i = counter;
408
409 if ((i<=(N-1)) )
410 {
411 ++counter;
412 return V[i];
413 }
414 else
415 {
416 sumtot = iterate_raw_vec(V, sumtot);
417 counter=2;
418 return V[1];
419 }
420 }
421
422 myuint_t MixMaxRng::precalc()
423 {
424 myuint_t temp = 0;
425 for (int i=0; i < N; ++i) { temp = MIXMAX_M
426 sumtot = temp;
427 return temp;
428 }
429
430 void MixMaxRng::state_init()
431 {
432 for (int i=1; i < N; ++i) { V[i] = 0; }
433 V[0] = 1;
434 counter = N; // set the counter to N if it
435 sumtot = 1;
436 }
437
438 void MixMaxRng::seed_spbox(myuint_t seed)
439 {
440 // a 64-bit LCG from Knuth line 26, in comb
441
442 if (seed == 0)
443 throw std::runtime_error("try seeding wit
444
445 const myuint_t MULT64=6364136223846793005UL
446 sumtot = 0;
447
448 myuint_t l = seed;
449
450 for (int i=0; i < N; ++i)
451 {
452 l*=MULT64; l = (l << 32) ^ (l>>32);
453 V[i] = l & M61;
454 sumtot = MIXMAX_MOD_MERSENNE(sumtot + V[(
455 }
456 counter = N; // set the counter to N if it
457 }
458
459 void MixMaxRng::seed_uniquestream( myID_t clus
460 {
461 state_init();
462 sumtot = apply_bigskip(V, V, clusterID, mac
463 counter = 1;
464 }
465
466 myuint_t MixMaxRng::apply_bigskip( myuint_t* V
467 {
468 /*
469 makes a derived state vector, Vout, from t
470 by skipping a large number of steps, deter
471
472 it is mathematically guaranteed that the s
473 1) at least one bit of ID is different
474 2) less than 10^100 numbers are drawn from
475 (this is good enough : a single CPU will n
476 even if it had a clock cycle of Planch tim
477
478 Caution: never apply this to a derived vec
479 and use it in all your runs, just change r
480
481 clusterID and machineID are provided for t
482 which is running in parallel on a large nu
483
484 did i repeat it enough times? the non-coll
485
486 */
487
488 const myuint_t skipMat17[128][17] =
489 #include "CLHEP/Random/mixmax_skip_N17.icc"
490 ;
491
492 const myuint_t* skipMat[128];
493 for (int i=0; i<128; ++i) { skipMat[i] = sk
494
495 myID_t IDvec[4] = {streamID, runID, machine
496 int r,i,j, IDindex;
497 myID_t id;
498 myuint_t Y[N], cum[N];
499 myuint_t coeff;
500 myuint_t* rowPtr;
501 myuint_t insumtot=0;
502
503 for (i=0; i<N; ++i) { Y[i] = Vin[i]; insumt
504 for (IDindex=0; IDindex<4; ++IDindex)
505 { // go from lower order to higher order ID
506 id=IDvec[IDindex];
507 //printf("now doing ID at level %d, with
508 r = 0;
509 while (id)
510 {
511 if (id & 1)
512 {
513 rowPtr = (myuint_t*)skipMat[r + IDind
514 for (i=0; i<N; ++i){ cum[i] = 0; }
515 for (j=0; j<N; ++j)
516 { // j is lag, enumerates terms of th
517 // for zero lag Y is already given
518 coeff = rowPtr[j]; // same coeff fo
519 for (i =0; i<N; ++i){
520 cum[i] = fmodmulM61( cum[i], coe
521 }
522 insumtot = iterate_raw_vec(Y, insum
523 }
524 insumtot=0;
525 for (i=0; i<N; ++i){ Y[i] = cum[i]; i
526 }
527 id = (id >> 1); ++r; // bring up the r-
528 }
529 }
530 insumtot=0;
531 for (i=0; i<N; ++i){ Vout[i] = Y[i]; insumt
532 // returns sumtot, and copy the vector over
533 return insumtot;
534 }
535
536 #if defined(__x86_64__)
537 myuint_t MixMaxRng::mod128(__uint128_t s)
538 {
539 myuint_t s1;
540 s1 = ( ( ((myuint_t)s)&M61 ) + ( ((myuin
541 return MIXMAX_MOD_MERSENNE(s1);
542 }
543 myuint_t MixMaxRng::fmodmulM61(myuint_t cum,
544 {
545 __uint128_t temp;
546 temp = (__uint128_t)a*(__uint128_t)b + cum
547 return mod128(temp);
548 }
549 #else // on all other platforms, including 32-
550 myuint_t MixMaxRng::fmodmulM61(myuint_t cum,
551 {
552 const myuint_t MASK32=0xFFFFFFFFULL;
553 myuint_t o,ph,pl,ah,al;
554 o=(s)*a;
555 ph = ((s)>>32);
556 pl = (s) & MASK32;
557 ah = a>>32;
558 al = a & MASK32;
559 o = (o & M61) + ((ph*ah)<<3) + ((ah*pl+al*
560 o += cum;
561 o = (o & M61) + ((o>>61));
562 return o;
563 }
564 #endif
565
566 void MixMaxRng::print_state() const
567 {
568 std::cout << "mixmax state, file version 1.
569 std::cout << "N=" << rng_get_N() << "; V[N]
570 for (int j=0; (j< (rng_get_N()-1) ); ++j) {
571 std::cout << V[rng_get_N()-1];
572 std::cout << "}; ";
573 std::cout << "counter= " << counter;
574 std::cout << "sumtot= " << sumtot << "\n";
575 }
576
577 MixMaxRng MixMaxRng::Branch()
578 {
579 sumtot = iterate_raw_vec(V, sumtot); counte
580 MixMaxRng tmp=*this;
581 tmp.BranchInplace(0); // daughter id
582 return tmp;
583 }
584
585 void MixMaxRng::BranchInplace(int id)
586 {
587 // Dont forget to iterate the mother, when
588 // a 64-bit LCG from Knuth line 26, is used
589 constexpr myuint_t MULT64=63641362238467930
590 myuint_t tmp=V[id];
591 V[1] *= MULT64; V[id] &= M61;
592 sumtot = MIXMAX_MOD_MERSENNE( sumtot + V[id
593 sumtot = iterate_raw_vec(V, sumtot);
594 counter = 1;
595 }
596
597 } // namespace CLHEP
598