41 const unsigned long MASK32=0xffffffff;
55 int numEngines = ++numberOfEngines;
56 setSeed(static_cast<long>(numEngines));
88 if (
this == &rng) {
return *
this; }
92 HepRandomEngine::operator=(rng);
104 FILE *fh= fopen(filename,
"w");
108 fprintf(fh,
"mixmax state, file version 1.0\n" );
109 fprintf(fh,
"N=%u; V[N]={",
rng_get_N() );
111 fprintf(fh,
"%llu, ",
S.
V[j] );
115 fprintf(fh,
"counter=%u; ",
S.
counter );
116 fprintf(fh,
"sumtot=%llu;\n",
S.
sumtot );
125 if( ( fin = fopen(filename,
"r") ) )
135 fprintf(stderr,
"mixmax -> read_state: error reading file %s\n", filename);
136 throw std::runtime_error(
"Error in reading state file");
141 if (!fscanf(fin,
"%llu", &
S.
V[0]) )
143 fprintf(stderr,
"mixmax -> read_state: error reading file %s\n", filename);
144 throw std::runtime_error(
"Error in reading state file");
150 if (!fscanf(fin,
", %llu", &vecVal) )
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");
161 fprintf(stderr,
"mixmax -> read_state: Invalid state vector value= %llu"
162 " ( must be less than %llu ) "
163 " obtained from reading file %s\n"
169 if (!fscanf( fin,
"}; counter=%i; ", &counter))
171 fprintf(stderr,
"mixmax -> read_state: error reading counter from file %s\n", filename);
172 throw std::runtime_error(
"Error in reading state file");
180 fprintf(stderr,
"mixmax -> read_state: Invalid counter = %d"
181 " Must be 0 <= counter < %u\n" , counter,
rng_get_N());
183 throw std::runtime_error(
"Error in reading state counter");
187 if (!fscanf( fin,
"sumtot=%llu\n", &sumtot))
189 fprintf(stderr,
"mixmax -> read_state: error reading checksum from file %s\n", filename);
190 throw std::runtime_error(
"Error in reading state file");
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");
203 std::cout << std::endl;
204 std::cout <<
"------- MixMaxRng engine status -------" << std::endl;
206 std::cout <<
" Current state vector is:" << std::endl;
208 std::cout <<
"---------------------------------------" << std::endl;
224 unsigned long seed0, seed1= 0, seed2= 0, seed3= 0;
227 seed0=
static_cast<unsigned long>(Seeds[0]) &
MASK32;
228 seed1=
static_cast<unsigned long>(Seeds[1]) &
MASK32;
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; }
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;
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__)
281 __asm__ __volatile__(
"pxor %0, %0;"
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++;};
326 for (
int i=0; i<size; ++i) { vect[i] =
flat(); }
329 MixMaxRng::operator double()
334 MixMaxRng::operator float()
336 return float(
flat() );
339 MixMaxRng::operator
unsigned int()
341 return static_cast<unsigned int>(get_next());
348 char beginMarker[] =
"MixMaxRng-begin";
349 char endMarker[] =
"MixMaxRng-end";
351 int pr = os.precision(24);
352 os << beginMarker <<
" ";
355 os <<
S.
V[i] <<
"\n";
359 os << endMarker <<
"\n";
366 std::vector<unsigned long>
v;
367 v.push_back (engineIDulong<MixMaxRng>());
370 v.push_back(static_cast<unsigned long>(
S.
V[i] &
MASK32));
372 v.push_back(static_cast<unsigned long>(
S.
V[i] >> 32 ));
375 v.push_back(static_cast<unsigned long>(
S.
counter));
377 v.push_back(static_cast<unsigned long>(
S.
sumtot >> 32));
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;
401 return "MixMaxRng-begin";
408 for (
int i=0; i<rng_get_N(); ++i) is >>
S.
V[i];
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";
422 std::cerr <<
"\nMixMaxRng::getState(): "
423 <<
"vector read wrong value of counter from file!"
424 <<
"\nInput stream is probably mispositioned now.\n";
429 std::cerr <<
"\nMixMaxRng::getState(): "
430 <<
"checksum disagrees with value stored in file!"
431 <<
"\nInput stream is probably mispositioned now.\n";
439 if ((v[0] & 0xffffffffUL) != engineIDulong<MixMaxRng>()) {
441 "\nMixMaxRng::get(): vector has wrong ID word - state unchanged\n";
451 "\nMixMaxRng::getState(): vector has wrong length - state unchanged\n";
462 std::cerr <<
"\nMixMaxRng::getState(): vector has wrong checksum!"
463 <<
"\nInput vector is probably mispositioned now.\n";
483 std::cerr <<
"MIXMAX ERROR: " <<
"Disallowed value of parameter N\n";
500 Y[0] = ( tempV = sumtotOld);
503 for (i=1; (i<
N); i++)
506 tempP =
modadd(tempP, Y[i]);
509 sumtot += tempV;
if (sumtot < tempV) {ovflow++;}
537 for (i=0; i <
N; i++){
555 for (i=0; i <
N; i++){
572 const myuint_t MULT64=6364136223846793005ULL;
576 if (seed == 0)
throw std::runtime_error(
"try seeding with nonzero seed next time");
580 for (i=0; i <
N; i++){
581 l*=MULT64; l = (l << 32) ^ (l>>32);
583 sumtot +=
S.
V[(i)];
if (sumtot <
S.
V[(i)]) {ovflow++;}
619 #include "CLHEP/Random/mixmax_skip_N17.icc"
623 for (
int i=0; i<128; i++) { skipMat[i] = skipMat17[i];}
625 myID_t IDvec[4] = {streamID, runID, machineID, clusterID};
633 for (i=0; i<
N; i++) { Y[i] = Vin[i]; sumtot =
modadd( sumtot, Vin[i]); } ;
634 for (IDindex=0; IDindex<4; IDindex++)
644 for (i=0; i<
N; i++){ cum[i] = 0; }
649 for (i =0; i<
N; i++){
655 for (i=0; i<
N; i++){ Y[i] = cum[i]; sumtot =
modadd( sumtot, cum[i]); } ;
661 for (i=0; i<
N; i++){ Vout[i] = Y[i]; sumtot =
modadd( sumtot, Y[i]); }
666 #if defined(__x86_64__)
667 myuint_t MixMaxRng::mod128(__uint128_t
s)
676 temp = (__uint128_t)a*(__uint128_t)b + cum;
679 #else // on all other platforms, including 32-bit linux, PPC and PPC64, ARM and all Windows
689 o = (o &
M61) + ((ph*ah)<<3) + ((ah*pl+al*ph + ((al*pl)>>32))>>29) ;
691 o = (o &
M61) + ((o>>61));
698 #if defined(__x86_64__) && defined(__GNUC__) && (!defined(__ICC))
702 __asm__ (
"addq %2, %0; "
717 std::cout <<
"mixmax state, file version 1.0\n";
718 std::cout <<
"N=" <<
rng_get_N() <<
"; V[N]={";
720 std::cout <<
S.
V[j] <<
", ";
724 std::cout <<
"counter= " <<
S.
counter;
725 std::cout <<
"sumtot= " <<
S.
sumtot <<
"\n";
740 constexpr
myuint_t MULT64=6364136223846793005ULL;
742 S.
V[1] *= MULT64;
S.
V[id] &=
M61;