Geant4  10.02.p02
mixmax.cc
Go to the documentation of this file.
1 // $Id:$
2 // -*- C++ -*-
3 //
4 // -----------------------------------------------------------------------
5 // MixMax Matrix PseudoRandom Number Generator
6 // --- MixMax ---
7 // class header file
8 // -----------------------------------------------------------------------
9 //
10 //
11 // Created by Konstantin Savvidy on Sun Feb 22 2004.
12 // As of version 0.99 and later, the code is being released under
13 // GNU Lesser General Public License v3
14 //
15 // Generator described in
16 // N.Z.Akopov, G.K.Savvidy and N.G.Ter-Arutyunian, Matrix Generator of Pseudorandom Numbers,
17 // J.Comput.Phys. 97, 573 (1991);
18 // Preprint EPI-867(18)-86, Yerevan Jun.1986;
19 //
20 // and
21 //
22 // K.Savvidy
23 // The MIXMAX random number generator
24 // Comp. Phys. Commun. (2015)
25 // http://dx.doi.org/10.1016/j.cpc.2015.06.003
26 //
27 // -----------------------------------------------------------------------
28 
29 #include <stdio.h>
30 #include <stdlib.h>
31 #include <stdint.h>
32 
33 #define __MIXMAX_C // do NOT define it in your own program, just include mixmax.h
34 
35 #include "CLHEP/Random/mixmax.h"
36 
37 namespace CLHEP {
38 
39 int iterate(rng_state_t* X){
40  X->sumtot = iterate_raw_vec(X->V, X->sumtot);
41  return 0;
42 }
43 
44 #if (SPECIALMUL!=0)
45 inline uint64_t MULWU (uint64_t k){ return (( (k)<<(SPECIALMUL) & M61) | ( (k) >> (BITS-SPECIALMUL)) ) ;}
46 #elif (SPECIALMUL==0)
47 inline uint64_t MULWU (uint64_t){ return 0;}
48 #else
49 #error SPECIALMUL not undefined
50 #endif
51 
52 myuint iterate_raw_vec(myuint* Y, myuint sumtotOld){
53  // operates with a raw vector, uses known sum of elements of Y
54  int i;
55 #ifdef SPECIAL
56  myuint temp2 = Y[1];
57 #endif
58  myuint tempP, tempV;
59  Y[0] = ( tempV = sumtotOld);
60  myuint sumtot = Y[0], ovflow = 0; // will keep a running sum of all new elements (except Y[0])
61  tempP = 0; // will keep a partial sum of all old elements (except Y[0])
62  for (i=1; i<N; i++){
63 #if (SPECIALMUL!=0)
64  myuint tempPO = MULWU(tempP);
65  tempP = modadd(tempP,Y[i]);
66  tempV = MOD_MERSENNE(tempV + tempP + tempPO); // edge cases ?
67 #else
68  tempP = modadd(tempP , Y[i]);
69  tempV = modadd(tempV , tempP);
70 #endif
71  Y[i] = tempV;
72  sumtot += tempV; if (sumtot < tempV) {ovflow++;}
73  }
74 #ifdef SPECIAL
75  temp2 = MOD_MULSPEC(temp2);
76  Y[2] = modadd( Y[2] , temp2 );
77  sumtot += temp2; if (sumtot < temp2) {ovflow++;}
78 #endif
79  return MOD_MERSENNE(MOD_MERSENNE(sumtot) + (ovflow <<3 ));
80 }
81 
82 myuint get_next(rng_state_t* X) {
83  return GET_BY_MACRO(X);
84 }
85 
86 double get_next_float(rng_state_t* X){
87  return get_next_float_BY_MACRO(X);
88 }
89 
90 void fill_array(rng_state_t* X, unsigned int n, double *array)
91 {
92  // Return an array of n random numbers uniformly distributed in (0,1]
93  unsigned int i,j;
94  const int M=N-1;
95  for (i=0; i<(n/M); i++){
96  iterate_and_fill_array(X, array+i*M);
97  }
98  unsigned int rem=(n % M);
99  if (rem) {
100  iterate(X);
101  for (j=0; j< (rem); j++){
102  array[M*i+j] = (int64_t)X->V[j] * (double)(INV_MERSBASE);
103  }
104  X->counter = j; // needed to continue with single fetches from the exact spot, but if you only use fill_array to get numbers then it is not necessary
105  }else{
106  X->counter = N;
107  }
108 }
109 
110 void iterate_and_fill_array(rng_state_t* X, double *array){
111  myuint* Y=X->V;
112  int i;
113  myuint tempP, tempV;
114 #if (SPECIAL != 0)
115  myuint temp2 = Y[1];
116 #endif
117  Y[0] = (tempV = modadd(Y[0] , X->sumtot));
118  //array[0] = (double)tempV * (double)(INV_MERSBASE);
119  myuint sumtot = 0, ovflow = 0; // will keep a running sum of all new elements (except Y[0])
120  tempP = 0; // will keep a partial sum of all old elements (except Y[0])
121  for (i=1; i<N; i++){
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);
126  }
127 #if (SPECIAL != 0)
128  temp2 = MOD_MULSPEC(temp2);
129  Y[2] = modadd( Y[2] , temp2 );
130  sumtot += temp2; if (sumtot < temp2) {ovflow++;}
131 #endif
132  X->sumtot = MOD_MERSENNE(MOD_MERSENNE(sumtot) + (ovflow <<3 ));
133 }
134 
135 myuint modadd(myuint foo, myuint bar){
136 #if defined(__x86_64__) && defined(USE_INLINE_ASM)
137  myuint out;
138  /* Assembler trick suggested by Andrzej Gòˆrlich */
139  __asm__ ("addq %2, %0; "
140  "btrq $61, %0; "
141  "adcq $0, %0; "
142  :"=r"(out)
143  :"0"(foo), "r"(bar)
144  );
145  return out;
146 #else
147  return MOD_MERSENNE(foo+bar);
148 #endif
149 }
150 
151 rng_state_t* rng_alloc()
152 {
153 /* allocate the state */
154  rng_state_t *p = (rng_state_t*)malloc(sizeof(rng_state_t));
155  p->fh=NULL; // by default, set the output file handle to stdout
156  return p;
157 }
158 
159 int rng_free(rng_state_t* X) /* free the memory occupied by the state */
160 {
161  free(X);
162  return 0;
163 }
164 
165 rng_state_t* rng_copy(myuint *Y)
166 {
167  /* copy the vector stored at Y, and return pointer to the newly allocated and initialized state.
168  It is the user's responsibility to make sure that Y is properly allocated with rng_alloc,
169  then pass Y->V or it can also be an array -- such as myuint Y[N+1] and Y[1]...Y[N] have been set to legal values [0 .. MERSBASE-1]
170  Partial sums on this new state are recalculated, and counter set to zero, so that when get_next is called,
171  it will output the initial vector before any new numbers are produced, call iterate(X) if you want to advance right away */
172  rng_state_t* X = rng_alloc();
173  myuint sumtot=0,ovflow=0;
174  X->counter = 2;
175  int i;
176  for ( i=0; i < N; i++){
177  X->V[i] = Y[i];
178  sumtot += X->V[(i)]; if (sumtot < X->V[(i)]) {ovflow++;}
179 
180  }
181  X->sumtot = MOD_MERSENNE(MOD_MERSENNE(sumtot) + (ovflow <<3 ));
182  return X;
183 }
184 
185 void seed_vielbein(rng_state_t* X, unsigned int index)
186 {
187 int i;
188  if (index<N){
189  for (i=0; i < N; i++){
190  X->V[i] = 0;
191  }
192  X->V[index] = 1;
193  }else{
194  fprintf(stderr, "Out of bounds index, is not ( 0 <= index < N )\n"); exit(ARRAY_INDEX_OUT_OF_BOUNDS);
195  }
196  X->counter = N; // set the counter to N if iteration should happen right away
197  //precalc(X);
198  X->sumtot = 1; //(index ? 1:0);
199  if (X->fh==NULL){X->fh=stdout;}
200 }
201 
202 void seed_spbox(rng_state_t* X, myuint seed)
203 { // a 64-bit LCG from Knuth line 26, in combination with a bit swap is used to seed
204  const myuint MULT64=6364136223846793005ULL;
205  int i;
206  myuint sumtot=0,ovflow=0;
207  if (seed == 0){
208  fprintf(stderr, " try seeding with nonzero seed next time!\n");
209  exit(SEED_WAS_ZERO);
210  }
211 
212  myuint l = seed;
213 
214  //X->V[0] = l & MERSBASE;
215  if (X->fh==NULL){X->fh=stdout;} // if the filehandle is not yet set, make it 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++;}
220  }
221  X->counter = N; // set the counter to N if iteration should happen right away
222  X->sumtot = MOD_MERSENNE(MOD_MERSENNE(sumtot) + (ovflow <<3 ));
223 }
224 
225 myuint precalc(rng_state_t* X){
226  int i;
227  myuint temp;
228  temp = 0;
229  for (i=0; i < N; i++){
230  temp = MOD_MERSENNE(temp + X->V[i]);
231  }
232  X->sumtot = temp;
233  return temp;
234 }
235 
236 
237 int rng_get_N(void){return N;}
238 
239 #if defined(__x86_64__)
240 inline myuint mod128(__uint128_t s){
241  myuint s1;
242  s1 = ( ( ((myuint)s)&MERSBASE ) + ( ((myuint)(s>>64)) * 8 ) + ( ((myuint)s) >>BITS) );
243  return MOD_MERSENNE(s1);
244 }
245 
246 inline myuint fmodmulM61(myuint cum, myuint a, myuint b){
247  __uint128_t temp;
248  temp = (__uint128_t)a*(__uint128_t)b + cum;
249  return mod128(temp);
250 }
251 
252 #else // on all other platforms, including 32-bit linux, PPC and PPC64 and all Windows
253 #define MASK32 0xFFFFFFFFULL
254 
255 inline myuint fmodmulM61(myuint cum, myuint s, myuint a)
256 {
257  register myuint o,ph,pl,ah,al;
258  o=(s)*a;
259  ph = ((s)>>32);
260  pl = (s) & MASK32;
261  ah = a>>32;
262  al = a & MASK32;
263  o = (o & M61) + ((ph*ah)<<3) + ((ah*pl+al*ph + ((al*pl)>>32))>>29) ;
264  o += cum;
265  o = (o & M61) + ((o>>61));
266  return o;
267 }
268 #endif
269 
270 void print_state(rng_state_t* X){
271  int j;
272  fprintf(X->fh, "mixmax state, file version 1.0\n" );
273  fprintf(X->fh, "N=%u; V[N]={", rng_get_N() );
274  for (j=0; (j< (rng_get_N()-1) ); j++) {
275  fprintf(X->fh, "%llu, ", X->V[j] );
276  }
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 );
281 }
282 
283 void read_state(rng_state_t* X, const char filename[] ){
284  // a function for reading the state from a file, after J. Apostolakis
285  FILE* fin;
286  if( ( fin = fopen(filename, "r") ) ){
287  char l=0;
288  while ( l != '{' ) { // 0x7B = "{"
289  l=fgetc(fin); // proceed until hitting opening bracket
290  }
291  ungetc(' ', fin);
292  }else{
293  fprintf(stderr, "mixmax -> read_state: error reading file %s\n", filename);
294  exit(ERROR_READING_STATE_FILE);
295  }
296 
297  myuint vecVal;
298  //printf("mixmax -> read_state: starting to read state from file\n");
299  if (!fscanf(fin, "%llu", &X->V[0]) ) {fprintf(stderr, "mixmax -> read_state: error reading file %s\n", filename); exit(ERROR_READING_STATE_FILE);}
300  //printf("V[%d] = %llu\n",0, X->V[0]);
301  int i;
302  for( i = 1; i < rng_get_N(); i++){
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);}
304  //printf("V[%d] = %llu\n",i, vecVal);
305  if( vecVal <= MERSBASE ){
306  X->V[i] = vecVal;
307  }else{
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);
312 
313  }
314  }
315 
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);}
318  if( counter <= N ) {
319  X->counter= counter;
320  }else{
321  fprintf(stderr, "mixmax -> read_state: Invalid counter = %d"
322  " Must be 0 <= counter < %u\n" , counter, N);
323  print_state(X);
324  exit(ERROR_READING_STATE_COUNTER);
325  }
326  precalc(X);
327  myuint sumtot;
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);}
329 
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);
333  }
334 // else{fprintf(stderr, "mixmax -> read_state: checksum ok: %llu == %llu\n",X->sumtot, sumtot);}
335  fclose(fin);
336 }
337 
338 
339 #define FUSEDMODMULVEC \
340 { for (i =0; i<N; i++){ \
341 cum[i] = fmodmulM61( cum[i], coeff , Y[i] ) ; \
342 } }
343 
344 
345 #define SKIPISON 1
346 
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 ){
349  seed_vielbein(Xin,0);
350  Xin->sumtot = apply_bigskip(Xin->V, Xin->V, clusterID, machineID, runID, streamID );
351  if (Xin->fh==NULL){Xin->fh=stdout;} // if the filehandle is not yet set, make it stdout
352 }
353 
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] );
356 }
357 
358 myuint apply_bigskip(myuint* Vout, myuint* Vin, myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID ){
359  /*
360  makes a derived state vector, Vout, from the mother state vector Vin
361  by skipping a large number of steps, determined by the given seeding ID's
362 
363  it is mathematically guaranteed that the substreams derived in this way from the SAME (!!!) Vin will not collide provided
364  1) at least one bit of ID is different
365  2) less than 10^100 numbers are drawn from the stream
366  (this is good enough : a single CPU will not exceed this in the lifetime of the universe, 10^19 sec,
367  even if it had a clock cycle of Planch time, 10^44 Hz )
368 
369  Caution: never apply this to a derived vector, just choose some mother vector Vin, for example the unit vector by seed_vielbein(X,0),
370  and use it in all your runs, just change runID to get completely nonoverlapping streams of random numbers on a different day.
371 
372  clusterID and machineID are provided for the benefit of large organizations who wish to ensure that a simulation
373  which is running in parallel on a large number of clusters and machines will have non-colliding source of random numbers.
374 
375  did i repeat it enough times? the non-collision guarantee is absolute, not probabilistic
376 
377  */
378 
379 
380  const myuint skipMat[128][N] =
381 
382 //#if (N==88)
383 //#include "CLHEP/Random/mixmax_skip_N88.icc" // to make this file, delete all except some chosen 128 rows of the coefficients table
384 //#elif (N==256)
386 //#elif (N==1000)
387 //#include "CLHEP/Random/mixmax_skip_N1000.icc"
388 //#elif (N==3150)
389 //#include "CLHEP/Random/mixmax_skip_N3150.icc"
390 //#endif
391  ;
392 
393  myID_t IDvec[4] = {streamID, runID, machineID, clusterID};
394  int r,i,j, IDindex;
395  myID_t id;
396  myuint Y[N], cum[N];
397  myuint coeff;
398  myuint* rowPtr;
399  myuint sumtot=0;
400 
401 
402  for (i=0; i<N; i++) { Y[i] = Vin[i]; sumtot = modadd( sumtot, Vin[i]); } ;
403  for (IDindex=0; IDindex<4; IDindex++) { // go from lower order to higher order ID
404  id=IDvec[IDindex];
405  //printf("now doing ID at level %d, with ID = %d\n", IDindex, id);
406  r = 0;
407  while (id){
408  if (id & 1) {
409  rowPtr = (myuint*)skipMat[r + IDindex*8*sizeof(myID_t)];
410  //printf("free coeff for row %d is %llu\n", r, rowPtr[0]);
411  for (i=0; i<N; i++){ cum[i] = 0; }
412  for (j=0; j<N; j++){ // j is lag, enumerates terms of the poly
413  // for zero lag Y is already given
414  coeff = rowPtr[j]; // same coeff for all i
416  sumtot = iterate_raw_vec(Y, sumtot);
417  }
418  sumtot=0;
419  for (i=0; i<N; i++){ Y[i] = cum[i]; sumtot = modadd( sumtot, cum[i]); } ;
420  }
421  id = (id >> 1); r++; // bring up the r-th bit in the ID
422  }
423  }
424  sumtot=0;
425  for (i=0; i<N; i++){ Vout[i] = Y[i]; sumtot = modadd( sumtot, Y[i]); } ; // returns sumtot, and copy the vector over to Vout
426  return (sumtot) ;
427 }
428 #else
429 #warning For this N, we dont have the skipping coefficients yet, using alternative method to seed
430 
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;
440  precalc(Xin);
441  Xin->sumtot = iterate_raw_vec(Xin->V, Xin->sumtot);
442  Xin->sumtot = iterate_raw_vec(Xin->V, Xin->sumtot);
443 }
444 #endif // SKIPISON
445 } // namespace CLHEP
446 
447 
448 
void print_state(rng_state_t *X)
Definition: mixmax.cc:270
double Y(double density)
int rng_get_N(void)
Definition: mixmax.cc:237
rng_state_t * rng_copy(myuint *Y)
Definition: mixmax.cc:165
myuint get_next(rng_state_t *X)
Definition: mixmax.cc:82
int rng_free(rng_state_t *X)
Definition: mixmax.cc:159
myuint fmodmulM61(myuint cum, myuint s, myuint a)
Definition: mixmax.cc:255
G4double a
Definition: TRTMaterials.hh:39
int iterate(rng_state_t *X)
Definition: mixmax.cc:39
static const double s
Definition: G4SIunits.hh:168
uint64_t MULWU(uint64_t)
Definition: mixmax.cc:47
void seed_spbox(rng_state_t *X, myuint seed)
Definition: mixmax.cc:202
double get_next_float(rng_state_t *X)
Definition: mixmax.cc:86
void seed_uniquestream(rng_state_t *Xin, myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID)
Definition: mixmax.cc:431
static const double bar
Definition: G4SIunits.hh:233
const G4int n
#define FUSEDMODMULVEC
Definition: mixmax.cc:339
#define BITS(n)
Definition: infback.cc:184
myuint iterate_raw_vec(myuint *Y, myuint sumtotOld)
Definition: mixmax.cc:52
void iterate_and_fill_array(rng_state_t *X, double *array)
Definition: mixmax.cc:110
rng_state_t * rng_alloc()
Definition: mixmax.cc:151
myuint precalc(rng_state_t *X)
Definition: mixmax.cc:225
myuint modadd(myuint foo, myuint bar)
Definition: mixmax.cc:135
void fill_array(rng_state_t *X, unsigned int n, double *array)
Definition: mixmax.cc:90
void read_state(rng_state_t *X, const char filename[])
Definition: mixmax.cc:283
void seed_vielbein(rng_state_t *X, unsigned int index)
Definition: mixmax.cc:185
#define MASK32
Definition: mixmax.cc:253