33 #define __MIXMAX_C // do NOT define it in your own program, just include mixmax.h
35 #include "CLHEP/Random/mixmax.h"
45 inline uint64_t
MULWU (uint64_t k){
return (( (k)<<(SPECIALMUL) & M61) | ( (k) >> (
BITS-SPECIALMUL)) ) ;}
47 inline uint64_t
MULWU (uint64_t){
return 0;}
49 #error SPECIALMUL not undefined
59 Y[0] = ( tempV = sumtotOld);
60 myuint sumtot = Y[0], ovflow = 0;
64 myuint tempPO =
MULWU(tempP);
65 tempP =
modadd(tempP,Y[i]);
66 tempV = MOD_MERSENNE(tempV + tempP + tempPO);
68 tempP =
modadd(tempP , Y[i]);
69 tempV =
modadd(tempV , tempP);
72 sumtot += tempV;
if (sumtot < tempV) {ovflow++;}
75 temp2 = MOD_MULSPEC(temp2);
76 Y[2] =
modadd( Y[2] , temp2 );
77 sumtot += temp2;
if (sumtot < temp2) {ovflow++;}
79 return MOD_MERSENNE(MOD_MERSENNE(sumtot) + (ovflow <<3 ));
83 return GET_BY_MACRO(X);
87 return get_next_float_BY_MACRO(X);
90 void fill_array(rng_state_t* X,
unsigned int n,
double *array)
95 for (i=0; i<(n/M); i++){
98 unsigned int rem=(n % M);
101 for (j=0; j< (rem); j++){
102 array[M*i+j] = (int64_t)X->V[j] * (
double)(INV_MERSBASE);
117 Y[0] = (tempV =
modadd(Y[0] , X->sumtot));
119 myuint sumtot = 0, ovflow = 0;
122 tempP =
modadd(tempP,Y[i]);
123 Y[i] = ( tempV =
modadd(tempV,tempP) );
124 sumtot += tempV;
if (sumtot < tempV) {ovflow++;}
125 array[i-1] = (int64_t)tempV * (
double)(INV_MERSBASE);
128 temp2 = MOD_MULSPEC(temp2);
129 Y[2] =
modadd( Y[2] , temp2 );
130 sumtot += temp2;
if (sumtot < temp2) {ovflow++;}
132 X->sumtot = MOD_MERSENNE(MOD_MERSENNE(sumtot) + (ovflow <<3 ));
136 #if defined(__x86_64__) && defined(USE_INLINE_ASM)
139 __asm__ (
"addq %2, %0; "
147 return MOD_MERSENNE(foo+bar);
154 rng_state_t *p = (rng_state_t*)malloc(
sizeof(rng_state_t));
173 myuint sumtot=0,ovflow=0;
176 for ( i=0; i < N; i++){
178 sumtot += X->V[(i)];
if (sumtot < X->V[(i)]) {ovflow++;}
181 X->sumtot = MOD_MERSENNE(MOD_MERSENNE(sumtot) + (ovflow <<3 ));
189 for (i=0; i < N; i++){
194 fprintf(stderr,
"Out of bounds index, is not ( 0 <= index < N )\n"); exit(ARRAY_INDEX_OUT_OF_BOUNDS);
199 if (X->fh==NULL){X->fh=stdout;}
204 const myuint MULT64=6364136223846793005ULL;
206 myuint sumtot=0,ovflow=0;
208 fprintf(stderr,
" try seeding with nonzero seed next time!\n");
215 if (X->fh==NULL){X->fh=stdout;}
216 for (i=0; i < N; i++){
217 l*=MULT64; l = (l << 32) ^ (l>>32);
218 X->V[i] = l & MERSBASE;
219 sumtot += X->V[(i)];
if (sumtot < X->V[(i)]) {ovflow++;}
222 X->sumtot = MOD_MERSENNE(MOD_MERSENNE(sumtot) + (ovflow <<3 ));
229 for (i=0; i < N; i++){
230 temp = MOD_MERSENNE(temp + X->V[i]);
239 #if defined(__x86_64__)
240 inline myuint mod128(__uint128_t
s){
242 s1 = ( ( ((myuint)s)&MERSBASE ) + ( ((myuint)(s>>64)) * 8 ) + ( ((myuint)
s) >>
BITS) );
243 return MOD_MERSENNE(s1);
246 inline myuint
fmodmulM61(myuint cum, myuint
a, myuint b){
248 temp = (__uint128_t)a*(__uint128_t)b + cum;
252 #else // on all other platforms, including 32-bit linux, PPC and PPC64 and all Windows
253 #define MASK32 0xFFFFFFFFULL
257 register myuint o,ph,pl,ah,al;
263 o = (o & M61) + ((ph*ah)<<3) + ((ah*pl+al*ph + ((al*pl)>>32))>>29) ;
265 o = (o & M61) + ((o>>61));
272 fprintf(X->fh,
"mixmax state, file version 1.0\n" );
273 fprintf(X->fh,
"N=%u; V[N]={",
rng_get_N() );
275 fprintf(X->fh,
"%llu, ", X->V[j] );
277 fprintf(X->fh,
"%llu", X->V[
rng_get_N()-1] );
278 fprintf(X->fh,
"}; " );
279 fprintf(X->fh,
"counter=%u; ", X->counter );
280 fprintf(X->fh,
"sumtot=%llu;\n", X->sumtot );
286 if( ( fin = fopen(filename,
"r") ) ){
293 fprintf(stderr,
"mixmax -> read_state: error reading file %s\n", filename);
294 exit(ERROR_READING_STATE_FILE);
299 if (!fscanf(fin,
"%llu", &X->V[0]) ) {fprintf(stderr,
"mixmax -> read_state: error reading file %s\n", filename); exit(ERROR_READING_STATE_FILE);}
303 if (!fscanf(fin,
", %llu", &vecVal) ) {fprintf(stderr,
"mixmax -> read_state: error reading vector component i=%d from file %s\n", i, filename); exit(ERROR_READING_STATE_FILE);}
305 if( vecVal <= MERSBASE ){
308 fprintf(stderr,
"mixmax -> read_state: Invalid state vector value= %llu"
309 " ( must be less than %llu ) "
310 " obtained from reading file %s\n"
311 , vecVal, MERSBASE, filename);
316 unsigned int counter;
317 if (!fscanf( fin,
"}; counter=%u; ", &counter)){fprintf(stderr,
"mixmax -> read_state: error reading counter from file %s\n", filename); exit(ERROR_READING_STATE_FILE);}
321 fprintf(stderr,
"mixmax -> read_state: Invalid counter = %d"
322 " Must be 0 <= counter < %u\n" , counter, N);
324 exit(ERROR_READING_STATE_COUNTER);
328 if (!fscanf( fin,
"sumtot=%llu\n", &sumtot)){fprintf(stderr,
"mixmax -> read_state: error reading checksum from file %s\n", filename); exit(ERROR_READING_STATE_FILE);}
330 if (X->sumtot != sumtot) {
331 fprintf(stderr,
"mixmax -> checksum error while reading state from file %s - corrupted?\n", filename);
332 exit(ERROR_READING_STATE_CHECKSUM);
339 #define FUSEDMODMULVEC \
340 { for (i =0; i<N; i++){ \
341 cum[i] = fmodmulM61( cum[i], coeff , Y[i] ) ; \
347 #if (BITS==61 && SKIPISON!=0)
348 void seed_uniquestream( rng_state_t* Xin, myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID ){
350 Xin->sumtot = apply_bigskip(Xin->V, Xin->V, clusterID, machineID, runID, streamID );
351 if (Xin->fh==NULL){Xin->fh=stdout;}
354 void branch_inplace( rng_state_t* Xin, myID_t* IDvec ){
355 Xin->sumtot = apply_bigskip(Xin->V, Xin->V, IDvec[3], IDvec[2], IDvec[1], IDvec[0] );
358 myuint apply_bigskip(myuint* Vout, myuint* Vin, myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID ){
380 const myuint skipMat[128][N] =
393 myID_t IDvec[4] = {streamID, runID, machineID, clusterID};
402 for (i=0; i<N; i++) { Y[i] = Vin[i]; sumtot =
modadd( sumtot, Vin[i]); } ;
403 for (IDindex=0; IDindex<4; IDindex++) {
409 rowPtr = (myuint*)skipMat[r + IDindex*8*
sizeof(myID_t)];
411 for (i=0; i<N; i++){ cum[i] = 0; }
419 for (i=0; i<N; i++){ Y[i] = cum[i]; sumtot =
modadd( sumtot, cum[i]); } ;
425 for (i=0; i<N; i++){ Vout[i] = Y[i]; sumtot =
modadd( sumtot, Y[i]); } ;
429 #warning For this N, we dont have the skipping coefficients yet, using alternative method to seed
431 void seed_uniquestream( rng_state_t* Xin, myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID ){
432 Xin->V[0] = (myuint)clusterID;
433 Xin->V[1] = (myuint)machineID;
434 Xin->V[2] = (myuint)runID;
435 Xin->V[3] = (myuint)streamID;
436 Xin->V[4] = (myuint)clusterID << 5;
437 Xin->V[5] = (myuint)machineID << 7;
438 Xin->V[6] = (myuint)runID << 11;
439 Xin->V[7] = (myuint)streamID << 13;
void print_state(rng_state_t *X)
rng_state_t * rng_copy(myuint *Y)
myuint get_next(rng_state_t *X)
int rng_free(rng_state_t *X)
myuint fmodmulM61(myuint cum, myuint s, myuint a)
int iterate(rng_state_t *X)
void seed_spbox(rng_state_t *X, myuint seed)
double get_next_float(rng_state_t *X)
void seed_uniquestream(rng_state_t *Xin, myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID)
myuint iterate_raw_vec(myuint *Y, myuint sumtotOld)
void iterate_and_fill_array(rng_state_t *X, double *array)
rng_state_t * rng_alloc()
myuint precalc(rng_state_t *X)
myuint modadd(myuint foo, myuint bar)
void fill_array(rng_state_t *X, unsigned int n, double *array)
void read_state(rng_state_t *X, const char filename[])
void seed_vielbein(rng_state_t *X, unsigned int index)