ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MixMaxRng.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file MixMaxRng.cc
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-2017
31 // =======================================================================
32 
33 #include "CLHEP/Random/Random.h"
34 #include "CLHEP/Random/MixMaxRng.h"
37 
38 #include <string.h> // for strcmp
39 #include <cmath>
40 
41 const unsigned long MASK32=0xffffffff;
42 
43 namespace CLHEP {
44 
45 namespace {
46  // Number of instances with automatic seed selection
47  CLHEP_ATOMIC_INT_TYPE numberOfEngines(0);
48 }
49 
50 static const int MarkerLen = 64; // Enough room to hold a begin or end marker.
51 
54 {
55  int numEngines = ++numberOfEngines;
56  setSeed(static_cast<long>(numEngines));
57 }
58 
61 {
62  theSeed=seed;
63  setSeed(seed);
64 }
65 
66 MixMaxRng::MixMaxRng(std::istream& is)
68 {
69  get(is);
70 }
71 
73 {
74 }
75 
77  : HepRandomEngine(rng)
78 {
79  S.V = rng.S.V;
80  S.sumtot= rng.S.sumtot;
81  S.counter= rng.S.counter;
82 }
83 
85 {
86  // Check assignment to self
87  //
88  if (this == &rng) { return *this; }
89 
90  // Copy base class data
91  //
92  HepRandomEngine::operator=(rng);
93 
94  S.V = rng.S.V;
95  S.sumtot= rng.S.sumtot;
96  S.counter= rng.S.counter;
97 
98  return *this;
99 }
100 
101 void MixMaxRng::saveStatus( const char filename[] ) const
102 {
103  // Create a C file-handle or an object that can accept the same output
104  FILE *fh= fopen(filename, "w");
105  if( fh )
106  {
107  int j;
108  fprintf(fh, "mixmax state, file version 1.0\n" );
109  fprintf(fh, "N=%u; V[N]={", rng_get_N() );
110  for (j=0; (j< (rng_get_N()-1) ); j++) {
111  fprintf(fh, "%llu, ", S.V[j] );
112  }
113  fprintf(fh, "%llu", S.V[rng_get_N()-1] );
114  fprintf(fh, "}; " );
115  fprintf(fh, "counter=%u; ", S.counter );
116  fprintf(fh, "sumtot=%llu;\n", S.sumtot );
117  fclose(fh);
118  }
119 }
120 
121 void MixMaxRng::restoreStatus( const char filename[] )
122 {
123  // a function for reading the state from a file
124  FILE* fin;
125  if( ( fin = fopen(filename, "r") ) )
126  {
127  char l=0;
128  while ( l != '{' ) { // 0x7B = "{"
129  l=fgetc(fin); // proceed until hitting opening bracket
130  }
131  ungetc(' ', fin);
132  }
133  else
134  {
135  fprintf(stderr, "mixmax -> read_state: error reading file %s\n", filename);
136  throw std::runtime_error("Error in reading state file");
137  }
138 
139  myuint_t vecVal;
140  //printf("mixmax -> read_state: starting to read state from file\n");
141  if (!fscanf(fin, "%llu", &S.V[0]) )
142  {
143  fprintf(stderr, "mixmax -> read_state: error reading file %s\n", filename);
144  throw std::runtime_error("Error in reading state file");
145  }
146 
147  int i;
148  for( i = 1; i < rng_get_N(); i++)
149  {
150  if (!fscanf(fin, ", %llu", &vecVal) )
151  {
152  fprintf(stderr, "mixmax -> read_state: error reading vector component i=%d from file %s\n", i, filename);
153  throw std::runtime_error("Error in reading state file");
154  }
155  if( vecVal <= MixMaxRng::M61 )
156  {
157  S.V[i] = vecVal;
158  }
159  else
160  {
161  fprintf(stderr, "mixmax -> read_state: Invalid state vector value= %llu"
162  " ( must be less than %llu ) "
163  " obtained from reading file %s\n"
164  , vecVal, MixMaxRng::M61, filename);
165  }
166  }
167 
168  int counter;
169  if (!fscanf( fin, "}; counter=%i; ", &counter))
170  {
171  fprintf(stderr, "mixmax -> read_state: error reading counter from file %s\n", filename);
172  throw std::runtime_error("Error in reading state file");
173  }
174  if( counter <= rng_get_N() )
175  {
176  S.counter= counter;
177  }
178  else
179  {
180  fprintf(stderr, "mixmax -> read_state: Invalid counter = %d"
181  " Must be 0 <= counter < %u\n" , counter, rng_get_N());
182  print_state();
183  throw std::runtime_error("Error in reading state counter");
184  }
185  precalc();
186  myuint_t sumtot;
187  if (!fscanf( fin, "sumtot=%llu\n", &sumtot))
188  {
189  fprintf(stderr, "mixmax -> read_state: error reading checksum from file %s\n", filename);
190  throw std::runtime_error("Error in reading state file");
191  }
192 
193  if (S.sumtot != sumtot)
194  {
195  fprintf(stderr, "mixmax -> checksum error while reading state from file %s - corrupted?\n", filename);
196  throw std::runtime_error("Error in reading state checksum");
197  }
198  fclose(fin);
199 }
200 
202 {
203  std::cout << std::endl;
204  std::cout << "------- MixMaxRng engine status -------" << std::endl;
205 
206  std::cout << " Current state vector is:" << std::endl;
207  print_state();
208  std::cout << "---------------------------------------" << std::endl;
209 }
210 
211 void MixMaxRng::setSeed(long longSeed, int /* extraSeed */)
212 {
213  //seed_uniquestream(0,0,0,longSeed);
214  theSeed = longSeed;
215  seed_spbox(longSeed);
216 }
217 
218 // Preferred Seeding method
219 // the values of 'Seeds' must be valid 32-bit integers
220 // Higher order bits will be ignored!!
221 
222 void MixMaxRng::setSeeds(const long* Seeds, int seedNum)
223 {
224  unsigned long seed0, seed1= 0, seed2= 0, seed3= 0;
225 
226  if( seedNum < 1 ) { // Assuming at least 2 seeds in vector...
227  seed0= static_cast<unsigned long>(Seeds[0]) & MASK32;
228  seed1= static_cast<unsigned long>(Seeds[1]) & MASK32;
229  }
230  else
231  {
232  if( seedNum < 4 ) {
233  seed0= static_cast<unsigned long>(Seeds[0]) & MASK32;
234  if( seedNum > 1){ seed1= static_cast<unsigned long>(Seeds[1]) & MASK32; }
235  if( seedNum > 2){ seed2= static_cast<unsigned long>(Seeds[2]) & MASK32; }
236  }
237  if( seedNum >= 4 ) {
238  seed0= static_cast<unsigned long>(Seeds[0]) & MASK32;
239  seed1= static_cast<unsigned long>(Seeds[1]) & MASK32;
240  seed2= static_cast<unsigned long>(Seeds[2]) & MASK32;
241  seed3= static_cast<unsigned long>(Seeds[3]) & MASK32;
242  }
243  }
244  theSeed = Seeds[0];
245  theSeeds = Seeds;
246  seed_uniquestream(seed3, seed2, seed1, seed0);
247 }
248 
250 {
251  return "MixMaxRng";
252 }
253 
254 constexpr int MixMaxRng::rng_get_N()
255 {
256  return N;
257 }
258 
259 constexpr long long int MixMaxRng::rng_get_SPECIAL()
260 {
261  return SPECIAL;
262 }
263 
265 {
266  return SPECIALMUL;
267 }
268 
269 double MixMaxRng::generate(int i)
270 {
271  S.counter++;
272 #if defined(__clang__) || defined(__llvm__)
273  return INV_M61*static_cast<double>(S.V[i]);
274 #elif defined(__GNUC__) && (__GNUC__ < 7) && (!defined(__ICC)) && defined(__x86_64__) && defined(__SSE2_MATH__)
275  int64_t Z=S.V[i];
276  double F=0.0;
277  //#warning Using the inline assembler
278  /* using SSE inline assemly to zero the xmm register, just before int64 -> double conversion,
279  not necessary in GCC-5 or better, but huge penalty on earlier compilers
280  */
281  __asm__ __volatile__( "pxor %0, %0;"
282  "cvtsi2sdq %1, %0;"
283  :"=x"(F)
284  :"r"(Z)
285  );
286  return F*INV_M61;
287 #else
288  //#warning other method
289  return convert1double(S.V[i]); //get_next_float_packbits();
290 #endif
291 }
292 
294 {
295  myuint_t* Y=S.V.data();
296  myuint_t tempP, tempV;
297  Y[0] = ( tempV = S.sumtot);
298  myuint_t sumtot = Y[0], ovflow = 0; // will keep a running sum of all new elements
299  tempP = 0; // will keep a partial sum of all old elements
300  myuint_t tempPO;
301  tempPO = MULWU(tempP); tempP = modadd(tempP, Y[1] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[1] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
302  tempPO = MULWU(tempP); tempP = modadd(tempP, Y[2] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[2] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
303  tempPO = MULWU(tempP); tempP = modadd(tempP, Y[3] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[3] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
304  tempPO = MULWU(tempP); tempP = modadd(tempP, Y[4] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[4] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
305  tempPO = MULWU(tempP); tempP = modadd(tempP, Y[5] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[5] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
306  tempPO = MULWU(tempP); tempP = modadd(tempP, Y[6] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[6] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
307  tempPO = MULWU(tempP); tempP = modadd(tempP, Y[7] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[7] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
308  tempPO = MULWU(tempP); tempP = modadd(tempP, Y[8] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[8] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
309  tempPO = MULWU(tempP); tempP = modadd(tempP, Y[9] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[9] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
310  tempPO = MULWU(tempP); tempP = modadd(tempP, Y[10]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[10] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
311  tempPO = MULWU(tempP); tempP = modadd(tempP, Y[11]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[11] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
312  tempPO = MULWU(tempP); tempP = modadd(tempP, Y[12]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[12] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
313  tempPO = MULWU(tempP); tempP = modadd(tempP, Y[13]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[13] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
314  tempPO = MULWU(tempP); tempP = modadd(tempP, Y[14]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[14] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
315  tempPO = MULWU(tempP); tempP = modadd(tempP, Y[15]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[15] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
316  tempPO = MULWU(tempP); tempP = modadd(tempP, Y[16]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[16] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
317  S.sumtot = MIXMAX_MOD_MERSENNE(MIXMAX_MOD_MERSENNE(sumtot) + (ovflow <<3 ));
318 
319  S.counter=2;
320  return double(S.V[1])*INV_M61;
321 }
322 
323 void MixMaxRng::flatArray(const int size, double* vect )
324 {
325  // fill_array( S, size, arrayDbl );
326  for (int i=0; i<size; ++i) { vect[i] = flat(); }
327 }
328 
329 MixMaxRng::operator double()
330 {
331  return flat();
332 }
333 
334 MixMaxRng::operator float()
335 {
336  return float( flat() );
337 }
338 
339 MixMaxRng::operator unsigned int()
340 {
341  return static_cast<unsigned int>(get_next());
342  // clhep_get_next returns a 64-bit integer, of which the lower 61 bits
343  // are random and upper 3 bits are zero
344 }
345 
346 std::ostream & MixMaxRng::put ( std::ostream& os ) const
347 {
348  char beginMarker[] = "MixMaxRng-begin";
349  char endMarker[] = "MixMaxRng-end";
350 
351  int pr = os.precision(24);
352  os << beginMarker << " ";
353  os << theSeed << "\n";
354  for (int i=0; i<rng_get_N(); ++i) {
355  os << S.V[i] << "\n";
356  }
357  os << S.counter << "\n";
358  os << S.sumtot << "\n";
359  os << endMarker << "\n";
360  os.precision(pr);
361  return os;
362 }
363 
364 std::vector<unsigned long> MixMaxRng::put () const
365 {
366  std::vector<unsigned long> v;
367  v.push_back (engineIDulong<MixMaxRng>());
368  for (int i=0; i<rng_get_N(); ++i)
369  {
370  v.push_back(static_cast<unsigned long>(S.V[i] & MASK32));
371  // little-ended order on all platforms
372  v.push_back(static_cast<unsigned long>(S.V[i] >> 32 ));
373  // pack uint64 into a data structure which is 32-bit on some platforms
374  }
375  v.push_back(static_cast<unsigned long>(S.counter));
376  v.push_back(static_cast<unsigned long>(S.sumtot & MASK32));
377  v.push_back(static_cast<unsigned long>(S.sumtot >> 32));
378  return v;
379 }
380 
381 std::istream & MixMaxRng::get ( std::istream& is)
382 {
383  char beginMarker [MarkerLen];
384  is >> std::ws;
385  is.width(MarkerLen); // causes the next read to the char* to be <=
386  // that many bytes, INCLUDING A TERMINATION \0
387  // (Stroustrup, section 21.3.2)
388  is >> beginMarker;
389  if (strcmp(beginMarker,"MixMaxRng-begin")) {
390  is.clear(std::ios::badbit | is.rdstate());
391  std::cerr << "\nInput stream mispositioned or"
392  << "\nMixMaxRng state description missing or"
393  << "\nwrong engine type found." << std::endl;
394  return is;
395  }
396  return getState(is);
397 }
398 
399 std::string MixMaxRng::beginTag ()
400 {
401  return "MixMaxRng-begin";
402 }
403 
404 std::istream & MixMaxRng::getState ( std::istream& is )
405 {
406  char endMarker[MarkerLen];
407  is >> theSeed;
408  for (int i=0; i<rng_get_N(); ++i) is >> S.V[i];
409  is >> S.counter;
410  myuint_t checksum;
411  is >> checksum;
412  is >> std::ws;
413  is.width(MarkerLen);
414  is >> endMarker;
415  if (strcmp(endMarker,"MixMaxRng-end")) {
416  is.clear(std::ios::badbit | is.rdstate());
417  std::cerr << "\nMixMaxRng state description incomplete."
418  << "\nInput stream is probably mispositioned now.\n";
419  return is;
420  }
421  if ( S.counter < 0 || S.counter > rng_get_N() ) {
422  std::cerr << "\nMixMaxRng::getState(): "
423  << "vector read wrong value of counter from file!"
424  << "\nInput stream is probably mispositioned now.\n";
425  return is;
426  }
427  precalc();
428  if ( checksum != S.sumtot) {
429  std::cerr << "\nMixMaxRng::getState(): "
430  << "checksum disagrees with value stored in file!"
431  << "\nInput stream is probably mispositioned now.\n";
432  return is;
433  }
434  return is;
435 }
436 
437 bool MixMaxRng::get (const std::vector<unsigned long> & v)
438 {
439  if ((v[0] & 0xffffffffUL) != engineIDulong<MixMaxRng>()) {
440  std::cerr <<
441  "\nMixMaxRng::get(): vector has wrong ID word - state unchanged\n";
442  return false;
443  }
444  return getState(v);
445 }
446 
447 bool MixMaxRng::getState (const std::vector<unsigned long> & v)
448 {
449  if (v.size() != VECTOR_STATE_SIZE ) {
450  std::cerr <<
451  "\nMixMaxRng::getState(): vector has wrong length - state unchanged\n";
452  return false;
453  }
454  for (int i=1; i<2*rng_get_N() ; i=i+2) {
455  S.V[i/2]= ( (v[i] & MASK32) | ( (myuint_t)(v[i+1]) << 32 ) );
456  // unpack from a data structure which is 32-bit on some platforms
457  }
458  S.counter = v[2*rng_get_N()+1];
459  precalc();
460  if ( ( (v[2*rng_get_N()+2] & MASK32)
461  | ( (myuint_t)(v[2*rng_get_N()+3]) << 32 ) ) != S.sumtot) {
462  std::cerr << "\nMixMaxRng::getState(): vector has wrong checksum!"
463  << "\nInput vector is probably mispositioned now.\n";
464  return false;
465  }
466  return true;
467 }
468 
470 {
471  switch (N)
472  {
473  case 17:
474  return 0;
475  break;
476  case 8:
477  return 0;
478  break;
479  case 240:
480  return fmodmulM61( 0, SPECIAL , (k) );
481  break;
482  default:
483  std::cerr << "MIXMAX ERROR: " << "Disallowed value of parameter N\n";
484  std::terminate();
485  break;
486  }
487 }
488 
490 {
491  return (( (k)<<(SPECIALMUL) & M61) ^ ( (k) >> (BITS-SPECIALMUL)) );
492 }
493 
495 {
496  // operates with a raw vector, uses known sum of elements of Y
497  int i;
498 
499  myuint_t tempP, tempV;
500  Y[0] = ( tempV = sumtotOld);
501  myuint_t sumtot = Y[0], ovflow = 0; // will keep a running sum of all new elements
502  tempP = 0; // will keep a partial sum of all old elements
503  for (i=1; (i<N); i++)
504  {
505  myuint_t tempPO = MULWU(tempP);
506  tempP = modadd(tempP, Y[i]);
507  tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); // new Y[i] = old Y[i] + old partial * m
508  Y[i] = tempV;
509  sumtot += tempV; if (sumtot < tempV) {ovflow++;}
510  }
511  return MIXMAX_MOD_MERSENNE(MIXMAX_MOD_MERSENNE(sumtot) + (ovflow <<3 ));
512 }
513 
515 {
516  int i;
517  i=S.counter;
518 
519  if ((i<=(N-1)) )
520  {
521  S.counter++;
522  return S.V[i];
523  }
524  else
525  {
526  S.sumtot = iterate_raw_vec(S.V.data(), S.sumtot);
527  S.counter=2;
528  return S.V[1];
529  }
530 }
531 
533 {
534  int i;
535  myuint_t temp;
536  temp = 0;
537  for (i=0; i < N; i++){
538  temp = MIXMAX_MOD_MERSENNE(temp + S.V[i]);
539  }
540  S.sumtot = temp;
541  return temp;
542 }
543 
545 {
546  myuint_t Z=get_next();
547  return convert1double(Z);
548 }
549 
550 void MixMaxRng::seed_vielbein(unsigned int index)
551 {
552  int i;
553  if (index<N)
554  {
555  for (i=0; i < N; i++){
556  S.V[i] = 0;
557  }
558  S.V[index] = 1;
559  }
560  else
561  {
562  std::terminate();
563  }
564  S.counter = N; // set the counter to N if iteration should happen right away
565  S.sumtot = 1;
566 }
567 
569 {
570  // a 64-bit LCG from Knuth line 26, in combination with a bit swap is used to seed
571 
572  const myuint_t MULT64=6364136223846793005ULL;
573  int i;
574 
575  myuint_t sumtot=0,ovflow=0;
576  if (seed == 0) throw std::runtime_error("try seeding with nonzero seed next time");
577 
578  myuint_t l = seed;
579 
580  for (i=0; i < N; i++){
581  l*=MULT64; l = (l << 32) ^ (l>>32);
582  S.V[i] = l & M61;
583  sumtot += S.V[(i)]; if (sumtot < S.V[(i)]) {ovflow++;}
584  }
585  S.counter = N; // set the counter to N if iteration should happen right away
586  S.sumtot = MIXMAX_MOD_MERSENNE(MIXMAX_MOD_MERSENNE(sumtot) + (ovflow <<3 ));
587 }
588 
589 void MixMaxRng::seed_uniquestream( myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID )
590 {
591  seed_vielbein(0);
592  S.sumtot = apply_bigskip(S.V.data(), S.V.data(), clusterID, machineID, runID, streamID );
593  S.counter = 1;
594 }
595 
596 myuint_t MixMaxRng::apply_bigskip( myuint_t* Vout, myuint_t* Vin, myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID )
597 {
598  /*
599  makes a derived state vector, Vout, from the mother state vector Vin
600  by skipping a large number of steps, determined by the given seeding ID's
601 
602  it is mathematically guaranteed that the substreams derived in this way from the SAME (!!!) Vin will not collide provided
603  1) at least one bit of ID is different
604  2) less than 10^100 numbers are drawn from the stream
605  (this is good enough : a single CPU will not exceed this in the lifetime of the universe, 10^19 sec,
606  even if it had a clock cycle of Planch time, 10^44 Hz )
607 
608  Caution: never apply this to a derived vector, just choose some mother vector Vin, for example the unit vector by seed_vielbein(X,0),
609  and use it in all your runs, just change runID to get completely nonoverlapping streams of random numbers on a different day.
610 
611  clusterID and machineID are provided for the benefit of large organizations who wish to ensure that a simulation
612  which is running in parallel on a large number of clusters and machines will have non-colliding source of random numbers.
613 
614  did i repeat it enough times? the non-collision guarantee is absolute, not probabilistic
615 
616  */
617 
618  const myuint_t skipMat17[128][17] =
619  #include "CLHEP/Random/mixmax_skip_N17.icc"
620  ;
621 
622  const myuint_t* skipMat[128];
623  for (int i=0; i<128; i++) { skipMat[i] = skipMat17[i];}
624 
625  myID_t IDvec[4] = {streamID, runID, machineID, clusterID};
626  int r,i,j, IDindex;
627  myID_t id;
628  myuint_t Y[N], cum[N];
629  myuint_t coeff;
630  myuint_t* rowPtr;
631  myuint_t sumtot=0;
632 
633  for (i=0; i<N; i++) { Y[i] = Vin[i]; sumtot = modadd( sumtot, Vin[i]); } ;
634  for (IDindex=0; IDindex<4; IDindex++)
635  { // go from lower order to higher order ID
636  id=IDvec[IDindex];
637  //printf("now doing ID at level %d, with ID = %d\n", IDindex, id);
638  r = 0;
639  while (id)
640  {
641  if (id & 1)
642  {
643  rowPtr = (myuint_t*)skipMat[r + IDindex*8*sizeof(myID_t)];
644  for (i=0; i<N; i++){ cum[i] = 0; }
645  for (j=0; j<N; j++)
646  { // j is lag, enumerates terms of the poly
647  // for zero lag Y is already given
648  coeff = rowPtr[j]; // same coeff for all i
649  for (i =0; i<N; i++){
650  cum[i] = fmodmulM61( cum[i], coeff , Y[i] ) ;
651  }
652  sumtot = iterate_raw_vec(Y, sumtot);
653  }
654  sumtot=0;
655  for (i=0; i<N; i++){ Y[i] = cum[i]; sumtot = modadd( sumtot, cum[i]); } ;
656  }
657  id = (id >> 1); r++; // bring up the r-th bit in the ID
658  }
659  }
660  sumtot=0;
661  for (i=0; i<N; i++){ Vout[i] = Y[i]; sumtot = modadd( sumtot, Y[i]); }
662  // returns sumtot, and copy the vector over to Vout
663  return (sumtot) ;
664 }
665 
666 #if defined(__x86_64__)
667  myuint_t MixMaxRng::mod128(__uint128_t s)
668  {
669  myuint_t s1;
670  s1 = ( ( ((myuint_t)s)&M61 ) + ( ((myuint_t)(s>>64)) * 8 ) + ( ((myuint_t)s) >>BITS) );
671  return MIXMAX_MOD_MERSENNE(s1);
672  }
674  {
675  __uint128_t temp;
676  temp = (__uint128_t)a*(__uint128_t)b + cum;
677  return mod128(temp);
678  }
679 #else // on all other platforms, including 32-bit linux, PPC and PPC64, ARM and all Windows
681  {
682  const myuint_t MASK32=0xFFFFFFFFULL;
683  myuint_t o,ph,pl,ah,al;
684  o=(s)*a;
685  ph = ((s)>>32);
686  pl = (s) & MASK32;
687  ah = a>>32;
688  al = a & MASK32;
689  o = (o & M61) + ((ph*ah)<<3) + ((ah*pl+al*ph + ((al*pl)>>32))>>29) ;
690  o += cum;
691  o = (o & M61) + ((o>>61));
692  return o;
693  }
694 #endif
695 
697 {
698 #if defined(__x86_64__) && defined(__GNUC__) && (!defined(__ICC))
699  //#warning Using assembler routine in modadd
700  myuint_t out;
701  /* Assembler trick suggested by Andrzej Gòˆrlich */
702  __asm__ ("addq %2, %0; "
703  "btrq $61, %0; "
704  "adcq $0, %0; "
705  :"=r"(out)
706  :"0"(foo), "r"(bar)
707  );
708  return out;
709 #else
710  return MIXMAX_MOD_MERSENNE(foo+bar);
711 #endif
712 }
713 
715 {
716  int j;
717  std::cout << "mixmax state, file version 1.0\n";
718  std::cout << "N=" << rng_get_N() << "; V[N]={";
719  for (j=0; (j< (rng_get_N()-1) ); j++) {
720  std::cout << S.V[j] << ", ";
721  }
722  std::cout << S.V[rng_get_N()-1];
723  std::cout << "}; ";
724  std::cout << "counter= " << S.counter;
725  std::cout << "sumtot= " << S.sumtot << "\n";
726 }
727 
729 {
730  S.sumtot = iterate_raw_vec(S.V.data(), S.sumtot); S.counter = 1;
731  MixMaxRng tmp=*this;
732  tmp.BranchInplace(0); // daughter id
733  return tmp;
734 }
735 
737 {
738  // Dont forget to iterate the mother, when branching the daughter, or else will have collisions!
739  // a 64-bit LCG from Knuth line 26, is used to mangle a vector component
740  constexpr myuint_t MULT64=6364136223846793005ULL;
741  myuint_t tmp=S.V[id];
742  S.V[1] *= MULT64; S.V[id] &= M61;
743  S.sumtot = MIXMAX_MOD_MERSENNE( S.sumtot + S.V[id] - tmp + M61);
744  S.sumtot = iterate_raw_vec(S.V.data(), S.sumtot);// printf("iterating!\n");
745  S.counter = 1;
746 }
747 
748 } // namespace CLHEP