Logo ROOT   6.19/01
Reference Guide
mixmax.h
Go to the documentation of this file.
1 /*
2  * mixmax.c
3  * A Pseudo-Random Number Generator
4  *
5  * Created by Konstantin Savvidy.
6  *
7  * The code is released under GNU Lesser General Public License v3
8  *
9  * G.K.Savvidy and N.G.Ter-Arutyunian,
10  * On the Monte Carlo simulation of physical systems,
11  * J.Comput.Phys. 97, 566 (1991);
12  * Preprint EPI-865-16-86, Yerevan, Jan. 1986
13  *
14  * K.Savvidy
15  * The MIXMAX random number generator
16  * Comp. Phys. Commun. 196 (2015), pp 161–165
17  * http://dx.doi.org/10.1016/j.cpc.2015.06.003
18  *
19  */
20 
21 #ifndef MIXMAX_H_
22 #define MIXMAX_H_
23 
24 //#define USE_INLINE_ASM
25 
26 //#ifdef __cplusplus
27 //extern "C" {
28 //#endif
29 
30 #ifndef ROOT_MM_N
31 #define N 240
32 /* The currently recommended generator is the three-parameter MIXMAX with
33 N=240, s=487013230256099140, m=2^51+1
34 
35  Other newly recommended N are N=8, 17 and 240,
36  as well as the ordinary MIXMAX with N=256 and s=487013230256099064
37 
38  Since the algorithm is linear in N, the cost per number is almost independent of N.
39  */
40 #else
41 #define N ROOT_MM_N
42 #endif
43 
44 #ifndef __LP64__
45 typedef uint64_t myuint;
46 //#warning but no problem, 'myuint' is 'uint64_t'
47 #else
48 typedef unsigned long long int myuint;
49 //#warning but no problem, 'myuint' is 'unsigned long long int'
50 #endif
51 
53 {
56  int counter;
57  FILE* fh;
58 };
59 
60 typedef struct rng_state_st rng_state_t; // C struct alias
61 
62 int rng_get_N(void); // get the N programmatically, useful for checking the value for which the library was compiled
63 
64 rng_state_t *rng_alloc(); /* allocate the state */
65 int rng_free(rng_state_t* X); /* free memory occupied by the state */
66 rng_state_t *rng_copy(myuint *Y); /* init from vector, takes the vector Y,
67  returns pointer to the newly allocated and initialized state */
68 void read_state(rng_state_t* X, const char filename[] );
69 void print_state(rng_state_t* X);
70  int iterate(rng_state_t* X);
71  myuint iterate_raw_vec(myuint* Y, myuint sumtotOld);
72 
73 
74 // FUNCTIONS FOR SEEDING
75 
76 typedef uint32_t myID_t;
77 
78 void seed_uniquestream(rng_state_t* X, myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID );
79 /*
80  best choice: will make a state vector from which you can get at least 10^100 numbers
81  guaranteed mathematically to be non-colliding with any other stream prepared from another set of 32bit IDs,
82  so long as it is different by at least one bit in at least one of the four IDs
83  -- useful if you are running a parallel simulation with many clusters, many CPUs each
84  */
85 
86 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
87 
88 void seed_vielbein(rng_state_t* X, unsigned int i); // seeds with the i-th unit vector, i = 0..N-1, for testing only
89 
90 
91 
92 // FUNCTIONS FOR GETTING RANDOM NUMBERS
93 
94 #ifdef __MIXMAX_C
95  myuint get_next(rng_state_t* X); // returns 64-bit int, which is between 1 and 2^61-1 inclusive
96  double get_next_float(rng_state_t* X); // returns double precision floating point number in (0,1]
97 #endif //__MIXMAX_C
98 
99 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)
100 
101 void iterate_and_fill_array(rng_state_t* X, double *array); // fills the array with N numbers
102 
104 /* needed if the state has been changed by something other than iterate, but no worries, seeding functions call this for you when necessary */
105 myuint apply_bigskip(myuint* Vout, myuint* Vin, myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID );
106 // applies a skip of some number of steps calculated from the four IDs
107 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
108 
109 
110 #define BITS 61
111 
112 /* magic with Mersenne Numbers */
113 
114 #define M61 2305843009213693951ULL
115 
116  myuint modadd(myuint foo, myuint bar);
117 //myuint modmulM61(myuint s, myuint a);
119 
120 #define MERSBASE M61 //xSUFF(M61)
121 #define MOD_PAYNE(k) ((((k)) & MERSBASE) + (((k)) >> BITS) ) // slightly faster than my old way, ok for addition
122 #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
123 #define MOD_MERSENNE(k) MOD_PAYNE(k)
124 
125 //#define INV_MERSBASE (0x1p-61)
126 #define INV_MERSBASE (0.4336808689942017736029811203479766845703E-18)
127 //const double INV_MERSBASE=(0.4336808689942017736029811203479766845703E-18); // gives "duplicate symbol" error
128 
129 // the charpoly is irreducible for the combinations of N and SPECIAL
130 
131 #if (N==256)
132 #define SPECIALMUL 0
133 #ifdef USE_MIXMAX_256_NEW
134 // for 1.1
135 #define SPECIAL 487013230256099064 // s=487013230256099064, m=1 -- good old MIXMAX
136 #define MOD_MULSPEC(k) fmodmulM61( 0, SPECIAL , (k) )
137 #else
138 // for 1.0
139 #define SPECIAL -1
140 #define MOD_MULSPEC(k) (MERSBASE - (k));
141 #endif
142 
143 #elif (N==8)
144 #define SPECIALMUL 53 // m=2^53+1
145 #define SPECIAL 0
146 
147 #elif (N==17)
148 #define SPECIALMUL 36 // m=2^36+1, other valid possibilities are m=2^13+1, m=2^19+1, m=2^24+1
149 #define SPECIAL 0
150 
151 #elif (N==40)
152 #define SPECIALMUL 42 // m=2^42+1
153 #define SPECIAL 0
154 
155 #elif (N==60)
156 #define SPECIALMUL 52 // m=2^52+1
157 #define SPECIAL 0
158 
159 #elif (N==96)
160 #define SPECIALMUL 55 // m=2^55+1
161 #define SPECIAL 0
162 
163 #elif (N==120)
164 #define SPECIALMUL 51 // m=2^51+1 and a SPECIAL=+1 (!!!)
165 #define SPECIAL 1
166 #define MOD_MULSPEC(k) (k)
167 
168 #elif (N==240)
169 #define SPECIALMUL 51 // m=2^51+1 and a SPECIAL=487013230256099140
170 #define SPECIAL 487013230256099140ULL
171 #define MOD_MULSPEC(k) fmodmulM61( 0, SPECIAL , (k) )
172 
173 #elif (N==44851)
174 #define SPECIALMUL 0
175 #define SPECIAL -3
176 #define MOD_MULSPEC(k) MOD_MERSENNE(3*(MERSBASE-(k)))
177 
178 
179 #else
180 #warning Not a verified N, you are on your own!
181 #define SPECIALMUL 58
182 #define SPECIAL 0
183 
184 #endif // list of interesting N for modulus M61 ends here
185 
186 
187 #ifndef __MIXMAX_C // c++ can put code into header files, why cant we? (with the inline declaration, should be safe from duplicate-symbol error)
188 
189 #define get_next(X) GET_BY_MACRO(X)
190 #define get_next_float(X) get_next_float_BY_MACRO(X)
191 
192 #endif // __MIXMAX_C
193 
194 
196  int i;
197  i=X->counter;
198 
199  if (i<=(N-1) ){
200  X->counter++;
201  return X->V[i];
202  }else{
203  X->sumtot = iterate_raw_vec(X->V, X->sumtot);
204  X->counter=2;
205  return X->V[1];
206  }
207 }
208 
210  /* cast to signed int trick suggested by Andrzej Görlich */
211  int64_t Z=(int64_t)GET_BY_MACRO(X);
212  double F;
213 #if defined(__GNUC__) && (__GNUC__ < 5) && (!defined(__ICC)) && defined(__x86_64__) && defined(__SSE2_MATH__) && defined(USE_INLINE_ASM)
214 //#warning Using the inline assembler
215 /* using SSE inline assemly to zero the xmm register, just before int64 -> double conversion,
216  not really necessary in GCC-5 or better, but huge penalty on earlier compilers
217  */
218  __asm__ __volatile__("pxor %0, %0; "
219  :"=x"(F)
220  );
221 #endif
222  F=Z;
223  return F*INV_MERSBASE;
224 }
225 
226 
227 // ERROR CODES - exit() is called with these
228 #define ARRAY_INDEX_OUT_OF_BOUNDS 0xFF01
229 #define SEED_WAS_ZERO 0xFF02
230 #define ERROR_READING_STATE_FILE 0xFF03
231 #define ERROR_READING_STATE_COUNTER 0xFF04
232 #define ERROR_READING_STATE_CHECKSUM 0xFF05
233 
234 // #ifdef __cplusplus
235 // }
236 // #endif
237 
238 //#define HOOKUP_GSL 1
239 
240 #ifdef HOOKUP_GSL // if you need to use mixmax through GSL, pass -DHOOKUP_GSL=1 to the compiler
241 #ifndef __MIXMAX_C
242 
243 #include <gsl/gsl_rng.h>
244 unsigned long gsl_get_next(void *vstate);
245 double gsl_get_next_float(void *vstate);
246 void seed_for_gsl(void *vstate, unsigned long seed);
247 
248 static const gsl_rng_type mixmax_type =
249 {"MIXMAX", /* name */
250  MERSBASE, /* RAND_MAX */
251  1, /* RAND_MIN */
252  sizeof (rng_state_t),
253  &seed_for_gsl,
254  &gsl_get_next,
255  &gsl_get_next_float
256 };
257 
258 unsigned long gsl_get_next(void *vstate) {
259  rng_state_t* X= (rng_state_t*)vstate;
260  return (unsigned long)get_next(X);
261 }
262 
263 double gsl_get_next_float(void *vstate) {
264  rng_state_t* X= (rng_state_t*)vstate;
265  return ( (double)get_next(X)) * INV_MERSBASE;
266 }
267 
268 void seed_for_gsl(void *vstate, unsigned long seed){
269  rng_state_t* X= (rng_state_t*)vstate;
270  seed_spbox(X,(myuint)seed);
271 }
272 
273 const gsl_rng_type *gsl_rng_mixmax = &mixmax_type;
274 
275 
276 #endif // HOOKUP_GSL
277 #endif // not inside __MIXMAX_C
278 #endif // closing MIXMAX_H_
int rng_free(rng_state_t *X)
Definition: mixmax.icc:162
myuint apply_bigskip(myuint *Vout, myuint *Vin, myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID)
Definition: mixmax.icc:398
int iterate(rng_state_t *X)
Definition: mixmax.icc:34
myuint V[N]
Definition: mixmax.h:54
myuint precalc(rng_state_t *X)
Definition: mixmax.icc:228
void fill_array(rng_state_t *X, unsigned int n, double *array)
Definition: mixmax.icc:85
void iterate_and_fill_array(rng_state_t *X, double *array)
Definition: mixmax.icc:105
myuint iterate_raw_vec(myuint *Y, myuint sumtotOld)
Definition: mixmax.icc:47
void seed_uniquestream(rng_state_t *X, myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID)
Definition: mixmax.icc:361
rng_state_t * rng_alloc()
Definition: mixmax.icc:154
#define get_next_float(X)
Definition: mixmax.h:190
myuint fmodmulM61(myuint cum, myuint s, myuint a)
Definition: mixmax.icc:258
void seed_vielbein(rng_state_t *X, unsigned int i)
Definition: mixmax.icc:188
void read_state(rng_state_t *X, const char filename[])
Definition: mixmax.icc:286
uint64_t myuint
Definition: mixmax.h:45
#define N
static constexpr double s
static const gsl_rng_type mixmax_type
#define F(x, y, z)
auto * a
Definition: textangle.C:12
FILE * fh
Definition: mixmax.h:57
myuint modadd(myuint foo, myuint bar)
Definition: mixmax.icc:137
uint32_t myID_t
Definition: mixmax.h:76
int counter
Definition: mixmax.h:56
void seed_spbox(rng_state_t *X, myuint seed)
Definition: mixmax.icc:205
myuint sumtot
Definition: mixmax.h:55
#define MERSBASE
myuint GET_BY_MACRO(rng_state_t *X)
Definition: mixmax.h:195
void print_state(rng_state_t *X)
Definition: mixmax.icc:273
#define INV_MERSBASE
#define get_next(X)
Definition: mixmax.h:189
struct rng_state_st rng_state_t
Definition: mixmax.h:60
static constexpr double bar
rng_state_t * rng_copy(myuint *Y)
Definition: mixmax.icc:168
void branch_inplace(rng_state_t *Xin, myID_t *ID)
Definition: mixmax.icc:394
double get_next_float_BY_MACRO(rng_state_t *X)
Definition: mixmax.h:209
const Int_t n
Definition: legend1.C:16
const gsl_rng_type * gsl_rng_mixmax
int rng_get_N(void)
Definition: mixmax.icc:240