44const unsigned long MASK32=0xffffffff;
59 setSeed(
static_cast<long>(numEngines));
91 if (
this == &rng) {
return *
this; }
95 HepRandomEngine::operator=(rng);
107 FILE *fh= fopen(filename,
"w");
111 fprintf(fh,
"mixmax state, file version 1.0\n" );
112 fprintf(fh,
"N=%u; V[N]={",
rng_get_N() );
114 fprintf(fh,
"%llu, ",
S.
V[j] );
118 fprintf(fh,
"counter=%u; ",
S.
counter );
119 fprintf(fh,
"sumtot=%llu;\n",
S.
sumtot );
128 if( ( fin = fopen(filename,
"r") ) )
138 fprintf(stderr,
"mixmax -> read_state: error reading file %s\n", filename);
139 throw std::runtime_error(
"Error in reading state file");
144 if (!fscanf(fin,
"%llu", &
S.
V[0]) )
146 fprintf(stderr,
"mixmax -> read_state: error reading file %s\n", filename);
147 throw std::runtime_error(
"Error in reading state file");
153 if (!fscanf(fin,
", %llu", &vecVal) )
155 fprintf(stderr,
"mixmax -> read_state: error reading vector component i=%d from file %s\n", i, filename);
156 throw std::runtime_error(
"Error in reading state file");
164 fprintf(stderr,
"mixmax -> read_state: Invalid state vector value= %llu"
165 " ( must be less than %llu ) "
166 " obtained from reading file %s\n"
172 if (!fscanf( fin,
"}; counter=%i; ", &counter))
174 fprintf(stderr,
"mixmax -> read_state: error reading counter from file %s\n", filename);
175 throw std::runtime_error(
"Error in reading state file");
183 fprintf(stderr,
"mixmax -> read_state: Invalid counter = %d"
184 " Must be 0 <= counter < %u\n" , counter,
rng_get_N());
186 throw std::runtime_error(
"Error in reading state counter");
190 if (!fscanf( fin,
"sumtot=%llu\n", &sumtot))
192 fprintf(stderr,
"mixmax -> read_state: error reading checksum from file %s\n", filename);
193 throw std::runtime_error(
"Error in reading state file");
198 fprintf(stderr,
"mixmax -> checksum error while reading state from file %s - corrupted?\n", filename);
199 throw std::runtime_error(
"Error in reading state checksum");
206 std::cout << std::endl;
207 std::cout <<
"------- MixMaxRng engine status -------" << std::endl;
209 std::cout <<
" Current state vector is:" << std::endl;
211 std::cout <<
"---------------------------------------" << std::endl;
227 myID_t seed0, seed1= 0, seed2= 0, seed3= 0;
237 if( seedNum > 1){ seed1=
static_cast<myID_t>(Seeds[1]) &
MASK32; }
238 if( seedNum > 2){ seed2=
static_cast<myID_t>(Seeds[2]) &
MASK32; }
275#if defined(__clang__) || defined(__llvm__)
276 return INV_M61*
static_cast<double>(
S.
V[i]);
277#elif defined(__GNUC__) && (__GNUC__ < 7) && (!defined(__ICC)) && defined(__x86_64__) && defined(__SSE2_MATH__)
284 __asm__ __volatile__(
"pxor %0, %0;"
304 tempPO =
MULWU(tempP); tempP =
modadd(tempP,
Y[1] ); tempV =
MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO);
Y[1] = tempV; sumtot += tempV;
if (sumtot < tempV) {ovflow++;};
305 tempPO =
MULWU(tempP); tempP =
modadd(tempP,
Y[2] ); tempV =
MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO);
Y[2] = tempV; sumtot += tempV;
if (sumtot < tempV) {ovflow++;};
306 tempPO =
MULWU(tempP); tempP =
modadd(tempP,
Y[3] ); tempV =
MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO);
Y[3] = tempV; sumtot += tempV;
if (sumtot < tempV) {ovflow++;};
307 tempPO =
MULWU(tempP); tempP =
modadd(tempP,
Y[4] ); tempV =
MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO);
Y[4] = tempV; sumtot += tempV;
if (sumtot < tempV) {ovflow++;};
308 tempPO =
MULWU(tempP); tempP =
modadd(tempP,
Y[5] ); tempV =
MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO);
Y[5] = tempV; sumtot += tempV;
if (sumtot < tempV) {ovflow++;};
309 tempPO =
MULWU(tempP); tempP =
modadd(tempP,
Y[6] ); tempV =
MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO);
Y[6] = tempV; sumtot += tempV;
if (sumtot < tempV) {ovflow++;};
310 tempPO =
MULWU(tempP); tempP =
modadd(tempP,
Y[7] ); tempV =
MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO);
Y[7] = tempV; sumtot += tempV;
if (sumtot < tempV) {ovflow++;};
311 tempPO =
MULWU(tempP); tempP =
modadd(tempP,
Y[8] ); tempV =
MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO);
Y[8] = tempV; sumtot += tempV;
if (sumtot < tempV) {ovflow++;};
312 tempPO =
MULWU(tempP); tempP =
modadd(tempP,
Y[9] ); tempV =
MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO);
Y[9] = tempV; sumtot += tempV;
if (sumtot < tempV) {ovflow++;};
313 tempPO =
MULWU(tempP); tempP =
modadd(tempP,
Y[10]); tempV =
MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO);
Y[10] = tempV; sumtot += tempV;
if (sumtot < tempV) {ovflow++;};
314 tempPO =
MULWU(tempP); tempP =
modadd(tempP,
Y[11]); tempV =
MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO);
Y[11] = tempV; sumtot += tempV;
if (sumtot < tempV) {ovflow++;};
315 tempPO =
MULWU(tempP); tempP =
modadd(tempP,
Y[12]); tempV =
MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO);
Y[12] = tempV; sumtot += tempV;
if (sumtot < tempV) {ovflow++;};
316 tempPO =
MULWU(tempP); tempP =
modadd(tempP,
Y[13]); tempV =
MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO);
Y[13] = tempV; sumtot += tempV;
if (sumtot < tempV) {ovflow++;};
317 tempPO =
MULWU(tempP); tempP =
modadd(tempP,
Y[14]); tempV =
MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO);
Y[14] = tempV; sumtot += tempV;
if (sumtot < tempV) {ovflow++;};
318 tempPO =
MULWU(tempP); tempP =
modadd(tempP,
Y[15]); tempV =
MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO);
Y[15] = tempV; sumtot += tempV;
if (sumtot < tempV) {ovflow++;};
319 tempPO =
MULWU(tempP); tempP =
modadd(tempP,
Y[16]); tempV =
MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO);
Y[16] = tempV; sumtot += tempV;
if (sumtot < tempV) {ovflow++;};
329 for (
int i=0; i<size; ++i) { vect[i] =
flat(); }
332MixMaxRng::operator double()
337MixMaxRng::operator float()
339 return float(
flat() );
342MixMaxRng::operator
unsigned int()
344 return static_cast<unsigned int>(get_next());
351 char beginMarker[] =
"MixMaxRng-begin";
352 char endMarker[] =
"MixMaxRng-end";
354 int pr = os.precision(24);
355 os << beginMarker <<
" ";
358 os <<
S.
V[i] <<
"\n";
362 os << endMarker <<
"\n";
369 std::vector<unsigned long> v;
370 v.push_back (engineIDulong<MixMaxRng>());
373 v.push_back(
static_cast<unsigned long>(
S.
V[i] &
MASK32));
375 v.push_back(
static_cast<unsigned long>(
S.
V[i] >> 32 ));
378 v.push_back(
static_cast<unsigned long>(
S.
counter));
380 v.push_back(
static_cast<unsigned long>(
S.
sumtot >> 32));
392 if (strcmp(beginMarker,
"MixMaxRng-begin")) {
393 is.clear(std::ios::badbit | is.rdstate());
394 std::cerr <<
"\nInput stream mispositioned or"
395 <<
"\nMixMaxRng state description missing or"
396 <<
"\nwrong engine type found." << std::endl;
404 return "MixMaxRng-begin";
418 if (strcmp(endMarker,
"MixMaxRng-end")) {
419 is.clear(std::ios::badbit | is.rdstate());
420 std::cerr <<
"\nMixMaxRng state description incomplete."
421 <<
"\nInput stream is probably mispositioned now.\n";
425 std::cerr <<
"\nMixMaxRng::getState(): "
426 <<
"vector read wrong value of counter from file!"
427 <<
"\nInput stream is probably mispositioned now.\n";
432 std::cerr <<
"\nMixMaxRng::getState(): "
433 <<
"checksum disagrees with value stored in file!"
434 <<
"\nInput stream is probably mispositioned now.\n";
442 if ((v[0] & 0xffffffffUL) != engineIDulong<MixMaxRng>()) {
444 "\nMixMaxRng::get(): vector has wrong ID word - state unchanged\n";
454 "\nMixMaxRng::getState(): vector has wrong length - state unchanged\n";
465 std::cerr <<
"\nMixMaxRng::getState(): vector has wrong checksum!"
466 <<
"\nInput vector is probably mispositioned now.\n";
486 std::cerr <<
"MIXMAX ERROR: " <<
"Disallowed value of parameter N\n";
503 Y[0] = ( tempV = sumtotOld);
506 for (i=1; (i<
N); i++)
512 sumtot += tempV;
if (sumtot < tempV) {ovflow++;}
540 for (i=0; i <
N; i++){
558 for (i=0; i <
N; i++){
575 const myuint_t MULT64=6364136223846793005ULL;
579 if (seed == 0)
throw std::runtime_error(
"try seeding with nonzero seed next time");
583 for (i=0; i <
N; i++){
584 l*=MULT64; l = (l << 32) ^ (l>>32);
586 sumtot +=
S.
V[(i)];
if (sumtot <
S.
V[(i)]) {ovflow++;}
626 for (
int i=0; i<128; i++) { skipMat[i] = skipMat17[i];}
628 myID_t IDvec[4] = {streamID, runID, machineID, clusterID};
636 for (i=0; i<
N; i++) {
Y[i] = Vin[i]; sumtot =
modadd( sumtot, Vin[i]); } ;
637 for (IDindex=0; IDindex<4; IDindex++)
647 for (i=0; i<
N; i++){ cum[i] = 0; }
652 for (i =0; i<
N; i++){
658 for (i=0; i<
N; i++){
Y[i] = cum[i]; sumtot =
modadd( sumtot, cum[i]); } ;
664 for (i=0; i<
N; i++){ Vout[i] =
Y[i]; sumtot =
modadd( sumtot,
Y[i]); }
669#if defined(__x86_64__)
670 myuint_t MixMaxRng::mod128(__uint128_t
s)
679 temp = (__uint128_t)a*(__uint128_t)b + cum;
692 o = (o &
M61) + ((ph*ah)<<3) + ((ah*pl+
al*ph + ((
al*pl)>>32))>>29) ;
694 o = (o &
M61) + ((o>>61));
701#if defined(__x86_64__) && defined(__GNUC__) && (!defined(__ICC))
705 __asm__ (
"addq %2, %0; "
720 std::cout <<
"mixmax state, file version 1.0\n";
721 std::cout <<
"N=" <<
rng_get_N() <<
"; V[N]={";
723 std::cout <<
S.
V[j] <<
", ";
727 std::cout <<
"counter= " <<
S.
counter;
728 std::cout <<
"sumtot= " <<
S.
sumtot <<
"\n";
743 constexpr myuint_t MULT64=6364136223846793005ULL;
745 S.
V[1] *= MULT64;
S.
V[id] &=
M61;
G4double Y(G4double density)
const unsigned long MASK32
#define MIXMAX_MOD_MERSENNE(k)
#define CLHEP_ATOMIC_INT_TYPE
void restoreStatus(const char filename[]="MixMaxRngState.conf")
static constexpr int BITS
static constexpr long long int rng_get_SPECIAL()
static constexpr long long int SPECIALMUL
void flatArray(const int size, double *vect)
static constexpr long long int SPECIAL
myuint_t MOD_MULSPEC(myuint_t k)
myuint_t iterate_raw_vec(myuint_t *Y, myuint_t sumtotOld)
void setSeed(long seed, int dum=0)
myuint_t fmodmulM61(myuint_t cum, myuint_t s, myuint_t a)
void seed_spbox(myuint_t)
myuint_t modadd(myuint_t foo, myuint_t bar)
static std::string beginTag()
void saveStatus(const char filename[]="MixMaxRngState.conf") const
static constexpr double INV_M61
static constexpr myuint_t M61
static constexpr int rng_get_SPECIALMUL()
myuint_t MULWU(myuint_t k)
virtual std::istream & get(std::istream &is)
MixMaxRng & operator=(const MixMaxRng &rng)
virtual std::istream & getState(std::istream &is)
std::vector< unsigned long > put() const
void BranchInplace(int id)
double get_next_float_packbits()
void setSeeds(const long *seeds, int seedNum=0)
static constexpr int rng_get_N()
double convert1double(myuint_t u)
static constexpr unsigned int VECTOR_STATE_SIZE
void seed_uniquestream(myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID)
myuint_t apply_bigskip(myuint_t *Vout, myuint_t *Vin, myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID)
static std::string engineName()
void seed_vielbein(unsigned int i)
CLHEP_ATOMIC_INT_TYPE numberOfEngines(0)
static const int MarkerLen
unsigned long long int myuint_t
static constexpr double bar
static constexpr double s
const G4double al
Mysterious coefficient that appears in the wavefunctions.
std::array< myuint_t, N > V