Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
mixmax.h
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 // The code is 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 #ifndef CLHEP_MIXMAX_H_
30 #define CLHEP_MIXMAX_H_ 1
31 
32 #include <stdio.h>
33 #include <stdint.h>
34 
35 #define USE_INLINE_ASM YES
36 
37 namespace CLHEP {
38 
39 #ifdef __cplusplus
40 extern "C" {
41 #endif
42 
43 const int N = 17;
44 /* The currently recommended N are 3150, 1260, 508, 256, 240, 88, 17, 8
45  Since the algorithm is linear in N, the cost per number is
46  almost independent of N.
47  */
48 
49 #ifndef __LP64__
50 typedef uint64_t myuint;
51 //#warning but no problem, 'myuint' is 'uint64_t'
52 #else
53 typedef unsigned long long int myuint;
54 //#warning but no problem, 'myuint' is 'unsigned long long int'
55 #endif
56 
58 {
59  myuint V[N];
60  myuint sumtot;
61  int counter;
62  FILE* fh;
63 };
64 
65 typedef struct rng_state_st rng_state_t; // C struct alias
66 
67 int rng_get_N(void); // get the N programmatically, useful for checking the value for which the library was compiled
68 
69 rng_state_t *rng_alloc(); /* allocate the state */
70 int rng_free(rng_state_t* X); /* free memory occupied by the state */
71 rng_state_t *rng_copy(myuint *Y); /* init from vector, takes the vector Y,
72  returns pointer to the newly allocated and initialized state */
73 void read_state(rng_state_t* X, const char filename[] );
74 void print_state(rng_state_t* X);
75  int iterate(rng_state_t* X);
76  myuint iterate_raw_vec(myuint* Y, myuint sumtotOld);
77 
78 
79 // FUNCTIONS FOR SEEDING
80 
81 typedef uint32_t myID_t;
82 
83 void seed_uniquestream(rng_state_t* X, myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID);
84 /*
85  best choice: will make a state vector from which you can get at least 10^100 numbers
86  guaranteed mathematically to be non-colliding with any other stream prepared from another set of 32bit IDs,
87  so long as it is different by at least one bit in at least one of the four IDs
88  -- useful if you are running a parallel simulation with many clusters, many CPUs each
89  */
90 
91 void seed_spbox(rng_state_t* X, myuint seed); // non-linear method, makes certified unique vectors, probability for streams to collide is < 1/10^4600
92 
93 void seed_vielbein(rng_state_t* X, unsigned int i); // seeds with the i-th unit vector, i = 0..N-1, for testing only
94 
95 
96 
97 // FUNCTIONS FOR GETTING RANDOM NUMBERS
98 
99 #ifdef __MIXMAX_C
100  myuint get_next(rng_state_t* X); // returns 64-bit int, which is between 1 and 2^61-1 inclusive
101  double get_next_float(rng_state_t* X); // returns double precision floating point number in (0,1]
102 #endif //__MIXMAX_C
103 
104 void fill_array(rng_state_t* X, unsigned int n, double *array); // fastest method: set n to a multiple of N (e.g. n=256)
105 
106 void iterate_and_fill_array(rng_state_t* X, double *array); // fills the array with N numbers
107 
108 myuint precalc(rng_state_t* X);
109 /* needed if the state has been changed by something other than iterate, but no worries, seeding functions call this for you when necessary */
110 myuint apply_bigskip(myuint* Vout, myuint* Vin, myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID );
111 // applies a skip of some number of steps calculated from the four IDs
112 void branch_inplace( rng_state_t* Xin, myID_t* ID ); // almost the same as apply_bigskip, but in-place and from a vector of IDs
113 
114 
115 #define BITS 61
116 
117 /* magic with Mersenne Numbers */
118 
119 #define M61 2305843009213693951ULL
120 
121  myuint modadd(myuint foo, myuint bar);
122  myuint modmulM61(myuint s, myuint a);
123  myuint fmodmulM61(myuint cum, myuint s, myuint a);
124 
125 #define MERSBASE M61 //xSUFF(M61)
126 #define MOD_PAYNE(k) ((((k)) & MERSBASE) + (((k)) >> BITS) ) // slightly faster than my old way, ok for addition
127 #define MOD_REM(k) ((k) % MERSBASE ) // latest Intel CPU is supposed to do this in one CPU cycle, but on my machines it seems to be 20% slower than the best tricks
128 #define MOD_MERSENNE(k) MOD_PAYNE(k)
129 
130 #define INV_MERSBASE (0.43368086899420177360298E-18L)
131 
132 
133 // the charpoly is irreducible for the combinations of N and SPECIAL and has maximal period for N=508, 256, half period for 1260, and 1/12 period for 3150
134 
135 // #if (N==256)
136 // #define SPECIALMUL 0
137 // #define SPECIAL 487013230256099064ULL // s=487013230256099064, m=1 -- good old MIXMAX
138 // #define MOD_MULSPEC(k) fmodmulM61( 0, SPECIAL , (k) );
139 
140 // #elif (N==17)
141 #define SPECIALMUL 36 // m=2^36+1
142 /*
143 #elif (N==8)
144 #define SPECIALMUL 53 // m=2^53+1
145 
146 #elif (N==40)
147 #define SPECIALMUL 42 // m=2^42+1
148 
149 #elif (N==96)
150 #define SPECIALMUL 55 // m=2^55+1
151 
152 #elif (N==64)
153 #define SPECIALMUL 55 // m=2^55 (!!!) and m=2^37+2
154 
155 #elif (N==120)
156 #define SPECIALMUL 51 // m=2^51+1 and a SPECIAL=+1 (!!!)
157 #define SPECIAL 1
158 #define MOD_MULSPEC(k) (k);
159 
160 #else
161 #warning Not a verified N, you are on your own!
162 #define SPECIALMUL 58
163 
164 #endif // list of interesting N for modulus M61 ends here
165 */
166 
167 #ifndef __MIXMAX_C // c++ can put code into header files, why cant we? (with the inline declaration, should be safe from duplicate-symbol error)
168 
169 #define get_next(X) GET_BY_MACRO(X)
170 #define get_next_float(X) get_next_float_BY_MACRO(X)
171 
172 #endif // __MIXMAX_C
173 
174 inline myuint GET_BY_MACRO(rng_state_t* X) {
175  int i;
176  i=X->counter;
177 
178  if (i<=(N-1) ){
179  X->counter++;
180  return X->V[i];
181  }else{
182  X->sumtot = iterate_raw_vec(X->V, X->sumtot);
183  X->counter=2;
184  return X->V[1];
185  }
186 }
187 
188 
190  int64_t Z=(int64_t)get_next(X);
191 #if defined(__x86_64__) && defined(__SSE__) && defined(__AVX__) && defined(USE_INLINE_ASM)
192  double F;
193  __asm__ __volatile__( "pxor %0, %0;"
194  //"cvtsi2sdq %1, %0;"
195  :"=x"(F)
196  //:"r"(Z)
197  );
198  F=Z;
199  return F*INV_MERSBASE;
200 #else
201  return Z*INV_MERSBASE;
202 #endif
203  }
204 
205 
206 
207 // ERROR CODES - exit() is called with these
208 #define ARRAY_INDEX_OUT_OF_BOUNDS 0xFF01
209 #define SEED_WAS_ZERO 0xFF02
210 #define ERROR_READING_STATE_FILE 0xFF03
211 #define ERROR_READING_STATE_COUNTER 0xFF04
212 #define ERROR_READING_STATE_CHECKSUM 0xFF05
213 
214 #ifdef __cplusplus
215 }
216 #endif
217 
218 //#define HOOKUP_GSL 1
219 
220 #ifdef HOOKUP_GSL // if you need to use mixmax through GSL, pass -DHOOKUP_GSL=1 to the compiler
221 
222 #include <gsl/gsl_rng.h>
223 unsigned long gsl_get_next(void *vstate);
224 double gsl_get_next_float(void *vstate);
225 void seed_for_gsl(void *vstate, unsigned long seed);
226 
227 static const gsl_rng_type mixmax_type =
228 {"MIXMAX", /* name */
229  MERSBASE, /* RAND_MAX */
230  1, /* RAND_MIN */
231  sizeof (rng_state_t),
232  &seed_for_gsl,
233  &gsl_get_next,
234  &gsl_get_next_float
235 };
236 
237 unsigned long gsl_get_next(void *vstate) {
238  rng_state_t* X= (rng_state_t*)vstate;
239  return (unsigned long)get_next(X);
240 }
241 
242 double gsl_get_next_float(void *vstate) {
243  rng_state_t* X= (rng_state_t*)vstate;
244  return ( (double)get_next(X)) * INV_MERSBASE;
245 }
246 
247 void seed_for_gsl(void *vstate, unsigned long seed){
248  rng_state_t* X= (rng_state_t*)vstate;
249  seed_spbox(X,(myuint)seed);
250 }
251 
252 const gsl_rng_type *gsl_rng_ran3 = &mixmax_type;
253 
254 
255 #endif // HOOKUP_GSL
256 
257 } // namespace CLHEP
258 
259 #endif // closing CLHEP_MIXMAX_H_
const int N
Definition: mixmax.h:43
void print_state(rng_state_t *X)
Definition: mixmax.cc:270
double Y(double density)
int rng_get_N(void)
Definition: mixmax.cc:237
myuint GET_BY_MACRO(rng_state_t *X)
Definition: mixmax.h:174
myuint get_next(rng_state_t *X)
Definition: mixmax.cc:82
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
void seed_vielbein(rng_state_t *X, unsigned int i)
Definition: mixmax.cc:185
uint32_t myID_t
Definition: mixmax.h:81
uint64_t myuint
Definition: mixmax.h:50
int rng_free(rng_state_t *X)
Definition: mixmax.cc:159
myuint fmodmulM61(myuint cum, myuint s, myuint a)
Definition: mixmax.cc:255
void seed_uniquestream(rng_state_t *X, myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID)
Definition: mixmax.cc:435
int iterate(rng_state_t *X)
Definition: mixmax.cc:39
const XML_Char * s
Definition: expat.h:262
myuint V[N]
Definition: mixmax.h:59
void seed_spbox(rng_state_t *X, myuint seed)
Definition: mixmax.cc:202
myuint modmulM61(myuint s, myuint a)
double get_next_float(rng_state_t *X)
Definition: mixmax.cc:86
long seed
Definition: chem4.cc:68
rng_state_t * rng_alloc()
Definition: mixmax.cc:151
rng_state_t * rng_copy(myuint *Y)
Definition: mixmax.cc:165
const G4int n
myuint apply_bigskip(myuint *Vout, myuint *Vin, myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID)
#define MERSBASE
Definition: mixmax.h:125
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
void branch_inplace(rng_state_t *Xin, myID_t *ID)
static constexpr double bar
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
#define INV_MERSBASE
Definition: mixmax.h:130
double get_next_float_BY_MACRO(rng_state_t *X)
Definition: mixmax.h:189
void read_state(rng_state_t *X, const char filename[])
Definition: mixmax.cc:283
struct rng_state_st rng_state_t
Definition: mixmax.h:65