| Geant4 Cross Reference |
1 //
2 // -*- C++ -*-
3 //
4 // -----------------------------------------------------------------------
5 // HEP Random
6 // --- MixMaxRng ---
7 // class implementation file
8 // -----------------------------------------------------------------------
9 //
10 // This file interfaces the MixMax PseudoRandom Number Generator
11 // proposed by:
12 //
13 // G.K.Savvidy and N.G.Ter-Arutyunian,
14 // On the Monte Carlo simulation of physical systems,
15 // J.Comput.Phys. 97, 566 (1991);
16 // Preprint EPI-865-16-86, Yerevan, Jan. 1986
17 // http://dx.doi.org/10.1016/0021-9991(91)90015-D
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.003
23 //
24 // K.Savvidy and G.Savvidy
25 // "Spectrum and Entropy of C-systems. MIXMAX random number generator"
26 // Chaos, Solitons & Fractals, Volume 91, (2016) pp. 33-38
27 // http://dx.doi.org/10.1016/j.chaos.2016.05.003
28 //
29 // =======================================================================
30 // Implementation by Konstantin Savvidy - Copyright 2004-2023
31 // July 2023 - Updated class structure upon suggestions from Marco Barbone
32 // September 2023 - fix (re-)initialization from Gabriele Cosmo
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 selection
53 CLHEP_ATOMIC_INT_TYPE numberOfEngines(0);
54 }
55
56 static const int MarkerLen = 64; // Enough room to hold a begin or end marker.
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 MixMaxRng& rng)
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 filename[] ) const
108 {
109 // Create a C file-handle or an object that can accept the same output
110 FILE *fh= fopen(filename, "w");
111 if( fh )
112 {
113 int j;
114 fprintf(fh, "mixmax state, file version 1.0\n" );
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 long)V[j] );
118 }
119 fprintf(fh, "%llu", (unsigned long long)V[rng_get_N()-1] );
120 fprintf(fh, "}; " );
121 fprintf(fh, "counter=%u; ", counter );
122 fprintf(fh, "sumtot=%llu;\n", (unsigned long long)sumtot );
123 fclose(fh);
124 }
125 }
126
127 void MixMaxRng::restoreStatus( const char filename[] )
128 {
129 // a function for reading the state from a file
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 opening bracket
136 }
137 ungetc(' ', fin);
138 }
139 else
140 {
141 fprintf(stderr, "mixmax -> read_state: error reading file %s\n", filename);
142 throw std::runtime_error("Error in reading state file");
143 }
144
145 myuint_t vecVal;
146 //printf("mixmax -> read_state: starting to read state from file\n");
147 if (!fscanf(fin, "%llu", (unsigned long long*) &V[0]) )
148 {
149 fprintf(stderr, "mixmax -> read_state: error reading file %s\n", filename);
150 throw std::runtime_error("Error in reading state file");
151 }
152
153 for (int i = 1; i < rng_get_N(); ++i)
154 {
155 if (!fscanf(fin, ", %llu", (unsigned long long*) &vecVal) )
156 {
157 fprintf(stderr, "mixmax -> read_state: error reading vector component i=%d from file %s\n", i, filename);
158 throw std::runtime_error("Error in reading state file");
159 }
160 if( vecVal <= MixMaxRng::M61 )
161 {
162 V[i] = vecVal;
163 }
164 else
165 {
166 fprintf(stderr, "mixmax -> read_state: Invalid state vector value= %llu"
167 " ( must be less than %llu ) "
168 " obtained from reading file %s\n"
169 , (unsigned long long)vecVal, (unsigned long long)MixMaxRng::M61, filename);
170 }
171 }
172
173 int incounter;
174 if (!fscanf( fin, "}; counter=%i; ", &incounter))
175 {
176 fprintf(stderr, "mixmax -> read_state: error reading counter from file %s\n", filename);
177 throw std::runtime_error("Error in reading state file");
178 }
179 if( incounter <= rng_get_N() )
180 {
181 counter = incounter;
182 }
183 else
184 {
185 fprintf(stderr, "mixmax -> read_state: Invalid counter = %d"
186 " Must be 0 <= counter < %u\n" , counter, rng_get_N());
187 print_state();
188 throw std::runtime_error("Error in reading state counter");
189 }
190 precalc();
191 myuint_t insumtot;
192 if (!fscanf( fin, "sumtot=%llu\n", (unsigned long long*) &insumtot))
193 {
194 fprintf(stderr, "mixmax -> read_state: error reading checksum from file %s\n", filename);
195 throw std::runtime_error("Error in reading state file");
196 }
197
198 if (sumtot != insumtot)
199 {
200 fprintf(stderr, "mixmax -> checksum error while reading state from file %s - corrupted?\n", filename);
201 throw std::runtime_error("Error in reading state checksum");
202 }
203 fclose(fin);
204 }
205
206 void MixMaxRng::showStatus() const
207 {
208 std::cout << std::endl;
209 std::cout << "------- MixMaxRng engine status -------" << std::endl;
210
211 std::cout << " Current state vector is:" << std::endl;
212 print_state();
213 std::cout << "---------------------------------------" << std::endl;
214 }
215
216 // Preferred Seeding method
217 // the values of 'Seeds' must be valid 32-bit integers
218 // Higher order bits will be ignored!!
219
220 void MixMaxRng::setSeeds(const long* Seeds, int seedNum)
221 {
222 myID_t seed0, seed1= 0, seed2= 0, seed3= 0;
223
224 if( seedNum < 1 ) { // Assuming at least 2 seeds in vector...
225 seed0= static_cast<myID_t>(Seeds[0]) & MASK32;
226 seed1= static_cast<myID_t>(Seeds[1]) & MASK32;
227 }
228 else
229 {
230 if( seedNum < 4 ) {
231 seed0= static_cast<myID_t>(Seeds[0]) & MASK32;
232 if( seedNum > 1){ seed1= static_cast<myID_t>(Seeds[1]) & MASK32; }
233 if( seedNum > 2){ seed2= static_cast<myID_t>(Seeds[2]) & MASK32; }
234 }
235 if( seedNum >= 4 ) {
236 seed0= static_cast<myID_t>(Seeds[0]) & MASK32;
237 seed1= static_cast<myID_t>(Seeds[1]) & MASK32;
238 seed2= static_cast<myID_t>(Seeds[2]) & MASK32;
239 seed3= static_cast<myID_t>(Seeds[3]) & MASK32;
240 }
241 }
242 theSeed = Seeds[0];
243 theSeeds = Seeds;
244 seed_uniquestream(seed3, seed2, seed1, seed0);
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, double* vect )
258 {
259 for (int i=0; i<size; ++i) { vect[i] = flat(); }
260 }
261
262 std::ostream & MixMaxRng::put ( std::ostream& os ) const
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 () const
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>(V[i] & MASK32));
287 // little-ended order on all platforms
288 vec.push_back(static_cast<unsigned long>(V[i] >> 32 ));
289 // pack uint64 into a data structure which is 32-bit on some platforms
290 }
291 vec.push_back(static_cast<unsigned long>(counter));
292 vec.push_back(static_cast<unsigned long>(sumtot & MASK32));
293 vec.push_back(static_cast<unsigned long>(sumtot >> 32));
294 return vec;
295 }
296
297 std::istream & MixMaxRng::get ( std::istream& is)
298 {
299 char beginMarker [MarkerLen];
300 is >> std::ws;
301 is.width(MarkerLen); // causes the next read to the char* to be <=
302 // that many bytes, INCLUDING A TERMINATION \0
303 // (Stroustrup, section 21.3.2)
304 is >> beginMarker;
305 if (strcmp(beginMarker,"MixMaxRng-begin")) {
306 is.clear(std::ios::badbit | is.rdstate());
307 std::cerr << "\nInput stream mispositioned or"
308 << "\nMixMaxRng state description missing or"
309 << "\nwrong engine type found." << std::endl;
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::istream& is )
321 {
322 char endMarker[MarkerLen];
323 is >> theSeed;
324 for (int i=0; i<rng_get_N(); ++i) is >> V[i];
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 description incomplete."
334 << "\nInput stream is probably mispositioned now.\n";
335 return is;
336 }
337 if ( counter < 0 || counter > rng_get_N() ) {
338 std::cerr << "\nMixMaxRng::getState(): "
339 << "vector read wrong value of counter from file!"
340 << "\nInput stream is probably mispositioned now.\n";
341 return is;
342 }
343 precalc();
344 if ( checksum != sumtot) {
345 std::cerr << "\nMixMaxRng::getState(): "
346 << "checksum disagrees with value stored in file!"
347 << "\nInput stream is probably mispositioned now.\n";
348 return is;
349 }
350 return is;
351 }
352
353 bool MixMaxRng::get (const std::vector<unsigned long> & vec)
354 {
355 if ((vec[0] & 0xffffffffUL) != engineIDulong<MixMaxRng>())
356 {
357 std::cerr <<
358 "\nMixMaxRng::get(): vector has wrong ID word - state unchanged\n";
359 return false;
360 }
361 return getState(vec);
362 }
363
364 bool MixMaxRng::getState (const std::vector<unsigned long> & vec)
365 {
366 if (vec.size() != VECTOR_STATE_SIZE ) {
367 std::cerr <<
368 "\nMixMaxRng::getState(): vector has wrong length - state unchanged\n";
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)(vec[i+1]) << 32 ) );
373 // unpack from a data structure which is 32-bit on some platforms
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]) << 32 ) ) != sumtot) {
379 std::cerr << "\nMixMaxRng::getState(): vector has wrong checksum!"
380 << "\nInput vector is probably mispositioned now.\n";
381 return false;
382 }
383 return true;
384 }
385
386 myuint_t MixMaxRng::iterate_raw_vec(myuint_t* Y, myuint_t sumtotOld)
387 {
388 // operates with a raw vector, uses known sum of elements of Y
389
390 myuint_t tempP, tempV;
391 Y[0] = ( tempV = sumtotOld);
392 myuint_t insumtot = Y[0], ovflow = 0; // will keep a running sum of all new elements
393 tempP = 0; // will keep a partial sum of all old elements
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+tempPO); // new Y[i] = old Y[i] + old partial * m
399 Y[i] = tempV;
400 insumtot += tempV; if (insumtot < tempV) {++ovflow;}
401 }
402 return MIXMAX_MOD_MERSENNE(MIXMAX_MOD_MERSENNE(insumtot) + (ovflow <<3 ));
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_MOD_MERSENNE(temp + V[i]); }
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 iteration should happen right away
435 sumtot = 1;
436 }
437
438 void MixMaxRng::seed_spbox(myuint_t seed)
439 {
440 // a 64-bit LCG from Knuth line 26, in combination with a bit swap is used to seed
441
442 if (seed == 0)
443 throw std::runtime_error("try seeding with nonzero seed next time");
444
445 const myuint_t MULT64=6364136223846793005ULL;
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[(i)]);
455 }
456 counter = N; // set the counter to N if iteration should happen right away
457 }
458
459 void MixMaxRng::seed_uniquestream( myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID )
460 {
461 state_init();
462 sumtot = apply_bigskip(V, V, clusterID, machineID, runID, streamID );
463 counter = 1;
464 }
465
466 myuint_t MixMaxRng::apply_bigskip( myuint_t* Vout, myuint_t* Vin, myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID )
467 {
468 /*
469 makes a derived state vector, Vout, from the mother state vector Vin
470 by skipping a large number of steps, determined by the given seeding ID's
471
472 it is mathematically guaranteed that the substreams derived in this way from the SAME (!!!) Vin will not collide provided
473 1) at least one bit of ID is different
474 2) less than 10^100 numbers are drawn from the stream
475 (this is good enough : a single CPU will not exceed this in the lifetime of the universe, 10^19 sec,
476 even if it had a clock cycle of Planch time, 10^44 Hz )
477
478 Caution: never apply this to a derived vector, just choose some mother vector Vin, for example the unit vector by state_init(),
479 and use it in all your runs, just change runID to get completely nonoverlapping streams of random numbers on a different day.
480
481 clusterID and machineID are provided for the benefit of large organizations who wish to ensure that a simulation
482 which is running in parallel on a large number of clusters and machines will have non-colliding source of random numbers.
483
484 did i repeat it enough times? the non-collision guarantee is absolute, not probabilistic
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] = skipMat17[i]; }
494
495 myID_t IDvec[4] = {streamID, runID, machineID, clusterID};
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]; insumtot = modadd( insumtot, Vin[i]); }
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 ID = %d\n", IDindex, id);
508 r = 0;
509 while (id)
510 {
511 if (id & 1)
512 {
513 rowPtr = (myuint_t*)skipMat[r + IDindex*8*sizeof(myID_t)];
514 for (i=0; i<N; ++i){ cum[i] = 0; }
515 for (j=0; j<N; ++j)
516 { // j is lag, enumerates terms of the poly
517 // for zero lag Y is already given
518 coeff = rowPtr[j]; // same coeff for all i
519 for (i =0; i<N; ++i){
520 cum[i] = fmodmulM61( cum[i], coeff , Y[i] ) ;
521 }
522 insumtot = iterate_raw_vec(Y, insumtot);
523 }
524 insumtot=0;
525 for (i=0; i<N; ++i){ Y[i] = cum[i]; insumtot = modadd( insumtot, cum[i]); } ;
526 }
527 id = (id >> 1); ++r; // bring up the r-th bit in the ID
528 }
529 }
530 insumtot=0;
531 for (i=0; i<N; ++i){ Vout[i] = Y[i]; insumtot = modadd( insumtot, Y[i]); }
532 // returns sumtot, and copy the vector over to Vout
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 ) + ( ((myuint_t)(s>>64)) * 8 ) + ( ((myuint_t)s) >>BITS) );
541 return MIXMAX_MOD_MERSENNE(s1);
542 }
543 myuint_t MixMaxRng::fmodmulM61(myuint_t cum, myuint_t a, myuint_t b)
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-bit linux, PPC and PPC64, ARM and all Windows
550 myuint_t MixMaxRng::fmodmulM61(myuint_t cum, myuint_t s, myuint_t a)
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*ph + ((al*pl)>>32))>>29) ;
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.0\n";
569 std::cout << "N=" << rng_get_N() << "; V[N]={";
570 for (int j=0; (j< (rng_get_N()-1) ); ++j) { std::cout << V[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); counter = 1;
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 branching the daughter, or else will have collisions!
588 // a 64-bit LCG from Knuth line 26, is used to mangle a vector component
589 constexpr myuint_t MULT64=6364136223846793005ULL;
590 myuint_t tmp=V[id];
591 V[1] *= MULT64; V[id] &= M61;
592 sumtot = MIXMAX_MOD_MERSENNE( sumtot + V[id] - tmp + M61);
593 sumtot = iterate_raw_vec(V, sumtot);
594 counter = 1;
595 }
596
597 } // namespace CLHEP
598