1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | |
8 | |
9 | |
10 | |
11 | |
12 | |
13 | |
14 | |
15 | |
16 | |
17 | |
18 | |
19 | |
20 | |
21 | |
22 | |
23 | |
24 | |
25 | |
26 | |
27 | |
28 | |
29 | |
30 | |
31 | |
32 | |
33 | |
34 | |
35 | |
36 | |
37 | #ifdef HAVE_CONFIG_H1 |
38 | #include <config.h> |
39 | #endif |
40 | |
41 | #include "random.h" |
42 | |
43 | #include <stdio.h> |
44 | #include <stdlib.h> |
45 | #ifdef HAVE_UNISTD_H |
46 | #include <unistd.h> |
47 | #endif |
48 | #include <time.h> |
49 | #include <math.h> |
50 | #ifdef GMX_NATIVE_WINDOWS |
51 | #include <process.h> |
52 | #endif |
53 | |
54 | #include "external/Random123-1.08/include/Random123/threefry.h" |
55 | |
56 | #include "gromacs/math/utilities.h" |
57 | #include "random_gausstable.h" |
58 | |
59 | #define RNG_N624 624 |
60 | #define RNG_M397 397 |
61 | #define RNG_MATRIX_A0x9908b0dfUL 0x9908b0dfUL /* constant vector a */ |
62 | #define RNG_UPPER_MASK0x80000000UL 0x80000000UL /* most significant w-r bits */ |
63 | #define RNG_LOWER_MASK0x7fffffffUL 0x7fffffffUL /* least significant r bits */ |
64 | |
65 | |
66 | |
67 | |
68 | |
69 | |
70 | |
71 | |
72 | #define GAUSS_TABLE14 14 /* the size of the gauss table is 2^GAUSS_TABLE */ |
73 | #define GAUSS_MASK((1 << 14) - 1) ((1 << GAUSS_TABLE14) - 1) |
74 | |
75 | |
76 | struct gmx_rng { |
77 | unsigned int mt[RNG_N624]; |
78 | int mti; |
79 | int has_saved; |
80 | double gauss_saved; |
81 | }; |
82 | |
83 | |
84 | |
85 | int |
86 | gmx_rng_n(void) |
87 | { |
88 | return RNG_N624; |
89 | } |
90 | |
91 | |
92 | gmx_rng_t |
93 | gmx_rng_init(unsigned int seed) |
94 | { |
95 | struct gmx_rng *rng; |
96 | |
97 | if ((rng = (struct gmx_rng *)malloc(sizeof(struct gmx_rng))) == NULL((void*)0)) |
98 | { |
99 | return NULL((void*)0); |
100 | } |
101 | |
102 | rng->has_saved = 0; |
103 | |
104 | rng->mt[0] = seed & 0xffffffffUL; |
105 | for (rng->mti = 1; rng->mti < RNG_N624; rng->mti++) |
106 | { |
107 | rng->mt[rng->mti] = |
108 | (1812433253UL * (rng->mt[rng->mti-1] ^ |
109 | (rng->mt[rng->mti-1] >> 30)) + rng->mti); |
110 | |
111 | |
112 | |
113 | |
114 | rng->mt[rng->mti] &= 0xffffffffUL; |
115 | |
116 | } |
117 | return rng; |
118 | } |
119 | |
120 | gmx_rng_t |
121 | gmx_rng_init_array(unsigned int seed[], int seed_length) |
122 | { |
123 | int i, j, k; |
124 | gmx_rng_t rng; |
125 | |
126 | if ((rng = gmx_rng_init(19650218UL)) == NULL((void*)0)) |
127 | { |
128 | return NULL((void*)0); |
129 | } |
130 | |
131 | i = 1; j = 0; |
132 | k = (RNG_N624 > seed_length ? RNG_N624 : seed_length); |
133 | for (; k; k--) |
134 | { |
135 | rng->mt[i] = (rng->mt[i] ^ ((rng->mt[i-1] ^ |
136 | (rng->mt[i-1] >> 30)) * 1664525UL)) |
137 | + seed[j] + j; |
138 | rng->mt[i] &= 0xffffffffUL; |
139 | i++; j++; |
140 | if (i >= RNG_N624) |
141 | { |
142 | rng->mt[0] = rng->mt[RNG_N624-1]; i = 1; |
143 | } |
144 | if (j >= seed_length) |
145 | { |
146 | j = 0; |
147 | } |
148 | } |
149 | for (k = RNG_N624-1; k; k--) |
150 | { |
151 | rng->mt[i] = (rng->mt[i] ^ ((rng->mt[i-1] ^ |
152 | (rng->mt[i-1] >> 30)) * |
153 | 1566083941UL)) |
154 | - i; |
155 | rng->mt[i] &= 0xffffffffUL; |
156 | i++; |
157 | if (i >= RNG_N624) |
158 | { |
159 | rng->mt[0] = rng->mt[RNG_N624-1]; i = 1; |
160 | } |
161 | } |
162 | |
163 | rng->mt[0] = 0x80000000UL; |
164 | |
165 | return rng; |
166 | } |
167 | |
168 | |
169 | void |
170 | gmx_rng_destroy(gmx_rng_t rng) |
171 | { |
172 | if (rng) |
173 | { |
174 | free(rng); |
175 | } |
176 | } |
177 | |
178 | |
179 | void |
180 | gmx_rng_get_state(gmx_rng_t rng, unsigned int *mt, int *mti) |
181 | { |
182 | int i; |
183 | |
184 | for (i = 0; i < RNG_N624; i++) |
185 | { |
186 | mt[i] = rng->mt[i]; |
187 | } |
188 | *mti = rng->mti; |
189 | } |
190 | |
191 | |
192 | void |
193 | gmx_rng_set_state(gmx_rng_t rng, unsigned int *mt, int mti) |
194 | { |
195 | int i; |
196 | |
197 | for (i = 0; i < RNG_N624; i++) |
198 | { |
199 | rng->mt[i] = mt[i]; |
200 | } |
201 | rng->mti = mti; |
202 | } |
203 | |
204 | |
205 | unsigned int |
206 | gmx_rng_make_seed(void) |
207 | { |
208 | FILE *fp; |
209 | unsigned int data; |
210 | int ret; |
211 | long my_pid; |
212 | |
213 | #ifdef HAVE_UNISTD_H |
214 | |
215 | |
216 | |
217 | |
218 | |
219 | |
220 | |
221 | |
222 | |
223 | |
224 | |
225 | |
226 | |
227 | fp = fopen("/dev/urandom", "rb"); |
228 | #else |
229 | fp = NULL((void*)0); |
230 | #endif |
231 | if (fp != NULL((void*)0)) |
232 | { |
233 | ret = fread(&data, sizeof(unsigned int), 1, fp); |
| Value stored to 'ret' is never read |
234 | fclose(fp); |
235 | } |
236 | else |
237 | { |
238 | |
239 | #ifdef GMX_NATIVE_WINDOWS |
240 | my_pid = (long)_getpid(); |
241 | #else |
242 | my_pid = (long)getpid(); |
243 | #endif |
244 | data = (unsigned int)(((long)time(NULL((void*)0))+my_pid) % (long)1000000); |
245 | } |
246 | return data; |
247 | } |
248 | |
249 | |
250 | |
251 | |
252 | |
253 | |
254 | static void |
255 | gmx_rng_update(gmx_rng_t rng) |
256 | { |
257 | unsigned int lastx, x1, x2, y, *mt; |
258 | int mti, k; |
259 | const unsigned int mag01[2] = {0x0UL, RNG_MATRIX_A0x9908b0dfUL}; |
260 | |
261 | |
262 | |
263 | mt = rng->mt; |
264 | mti = rng->mti; |
265 | |
266 | x1 = mt[0]; |
267 | for (k = 0; k < RNG_N624-RNG_M397-3; k += 4) |
268 | { |
269 | x2 = mt[k+1]; |
270 | y = (x1 & RNG_UPPER_MASK0x80000000UL) | (x2 & RNG_LOWER_MASK0x7fffffffUL); |
271 | mt[k] = mt[k+RNG_M397] ^ (y >> 1) ^ mag01[y & 0x1UL]; |
272 | x1 = mt[k+2]; |
273 | y = (x2 & RNG_UPPER_MASK0x80000000UL) | (x1 & RNG_LOWER_MASK0x7fffffffUL); |
274 | mt[k+1] = mt[k+RNG_M397+1] ^ (y >> 1) ^ mag01[y & 0x1UL]; |
275 | x2 = mt[k+3]; |
276 | y = (x1 & RNG_UPPER_MASK0x80000000UL) | (x2 & RNG_LOWER_MASK0x7fffffffUL); |
277 | mt[k+2] = mt[k+RNG_M397+2] ^ (y >> 1) ^ mag01[y & 0x1UL]; |
278 | x1 = mt[k+4]; |
279 | y = (x2 & RNG_UPPER_MASK0x80000000UL) | (x1 & RNG_LOWER_MASK0x7fffffffUL); |
280 | mt[k+3] = mt[k+RNG_M397+3] ^ (y >> 1) ^ mag01[y & 0x1UL]; |
281 | } |
282 | x2 = mt[k+1]; |
283 | y = (x1 & RNG_UPPER_MASK0x80000000UL) | (x2 & RNG_LOWER_MASK0x7fffffffUL); |
284 | mt[k] = mt[k+RNG_M397] ^ (y >> 1) ^ mag01[y & 0x1UL]; |
285 | k++; |
286 | x1 = mt[k+1]; |
287 | y = (x2 & RNG_UPPER_MASK0x80000000UL) | (x1 & RNG_LOWER_MASK0x7fffffffUL); |
288 | mt[k] = mt[k+RNG_M397] ^ (y >> 1) ^ mag01[y & 0x1UL]; |
289 | k++; |
290 | x2 = mt[k+1]; |
291 | y = (x1 & RNG_UPPER_MASK0x80000000UL) | (x2 & RNG_LOWER_MASK0x7fffffffUL); |
292 | mt[k] = mt[k+RNG_M397] ^ (y >> 1) ^ mag01[y & 0x1UL]; |
293 | k++; |
294 | for (; k < RNG_N624-1; k += 4) |
295 | { |
296 | x1 = mt[k+1]; |
297 | y = (x2 & RNG_UPPER_MASK0x80000000UL) | (x1 & RNG_LOWER_MASK0x7fffffffUL); |
298 | mt[k] = mt[k+(RNG_M397-RNG_N624)] ^ (y >> 1) ^ mag01[y & 0x1UL]; |
299 | x2 = mt[k+2]; |
300 | y = (x1 & RNG_UPPER_MASK0x80000000UL) | (x2 & RNG_LOWER_MASK0x7fffffffUL); |
301 | mt[k+1] = mt[k+(RNG_M397-RNG_N624)+1] ^ (y >> 1) ^ mag01[y & 0x1UL]; |
302 | x1 = mt[k+3]; |
303 | y = (x2 & RNG_UPPER_MASK0x80000000UL) | (x1 & RNG_LOWER_MASK0x7fffffffUL); |
304 | mt[k+2] = mt[k+(RNG_M397-RNG_N624)+2] ^ (y >> 1) ^ mag01[y & 0x1UL]; |
305 | x2 = mt[k+4]; |
306 | y = (x1 & RNG_UPPER_MASK0x80000000UL) | (x2 & RNG_LOWER_MASK0x7fffffffUL); |
307 | mt[k+3] = mt[k+(RNG_M397-RNG_N624)+3] ^ (y >> 1) ^ mag01[y & 0x1UL]; |
308 | } |
309 | y = (x2 & RNG_UPPER_MASK0x80000000UL) | (mt[0] & RNG_LOWER_MASK0x7fffffffUL); |
310 | mt[RNG_N624-1] = mt[RNG_M397-1] ^ (y >> 1) ^ mag01[y & 0x1UL]; |
311 | |
312 | rng->mti = 0; |
313 | } |
314 | |
315 | |
316 | real |
317 | gmx_rng_gaussian_real(gmx_rng_t rng) |
318 | { |
319 | real x, y, r; |
320 | |
321 | if (rng->has_saved) |
322 | { |
323 | rng->has_saved = 0; |
324 | return rng->gauss_saved; |
325 | } |
326 | else |
327 | { |
328 | do |
329 | { |
330 | x = 2.0*gmx_rng_uniform_real(rng)-1.0; |
331 | y = 2.0*gmx_rng_uniform_real(rng)-1.0; |
332 | r = x*x+y*y; |
333 | } |
334 | while (r > 1.0 || r == 0.0); |
335 | |
336 | r = sqrt(-2.0*log(r)/r); |
337 | rng->gauss_saved = y*r; |
338 | rng->has_saved = 1; |
339 | return x*r; |
340 | } |
341 | } |
342 | |
343 | |
344 | |
345 | |
346 | |
347 | |
348 | |
349 | |
350 | |
351 | unsigned int |
352 | gmx_rng_uniform_uint32(gmx_rng_t rng) |
353 | { |
354 | unsigned int y; |
355 | |
356 | if (rng->mti == RNG_N624) |
357 | { |
358 | gmx_rng_update(rng); |
359 | } |
360 | y = rng->mt[rng->mti++]; |
361 | |
362 | y ^= (y >> 11); |
363 | y ^= (y << 7) & 0x9d2c5680UL; |
364 | y ^= (y << 15) & 0xefc60000UL; |
365 | y ^= (y >> 18); |
366 | |
367 | return y; |
368 | } |
369 | |
370 | |
371 | |
372 | |
373 | |
374 | |
375 | real |
376 | gmx_rng_uniform_real(gmx_rng_t rng) |
377 | { |
378 | if (sizeof(real) == sizeof(double)) |
379 | { |
380 | return ((double)gmx_rng_uniform_uint32(rng))*(1.0/4294967296.0); |
381 | } |
382 | else |
383 | { |
384 | return ((float)gmx_rng_uniform_uint32(rng))*(1.0/4294967423.0); |
385 | } |
386 | |
387 | |
388 | |
389 | |
390 | |
391 | } |
392 | real |
393 | |
394 | gmx_rng_gaussian_table(gmx_rng_t rng) |
395 | { |
396 | unsigned int i; |
397 | |
398 | i = gmx_rng_uniform_uint32(rng); |
399 | |
400 | |
401 | return gaussian_table[i >> (32 - GAUSS_TABLE14)]; |
402 | } |
403 | |
404 | void |
405 | gmx_rng_cycle_2uniform(gmx_int64_t ctr1, gmx_int64_t ctr2, |
406 | gmx_int64_t key1, gmx_int64_t key2, |
407 | double* rnd) |
408 | { |
409 | const gmx_int64_t mask_53bits = 0x1FFFFFFFFFFFFF; |
410 | const double two_power_min53 = 1.0/9007199254740992.0; |
411 | |
412 | threefry2x64_ctr_t ctr = {{ctr1, ctr2}}; |
413 | threefry2x64_key_t key = {{key1, key2}}; |
414 | threefry2x64_ctr_t rand = threefry2x64(ctr, key)threefry2x64_R(threefry2x64_rounds, ctr, key); |
415 | |
416 | rnd[0] = (rand.v[0] & mask_53bits)*two_power_min53; |
417 | rnd[1] = (rand.v[1] & mask_53bits)*two_power_min53; |
418 | } |
419 | |
420 | void |
421 | gmx_rng_cycle_3gaussian_table(gmx_int64_t ctr1, gmx_int64_t ctr2, |
422 | gmx_int64_t key1, gmx_int64_t key2, |
423 | real* rnd) |
424 | { |
425 | threefry2x64_ctr_t ctr = {{ctr1, ctr2}}; |
426 | threefry2x64_key_t key = {{key1, key2}}; |
427 | threefry2x64_ctr_t rand = threefry2x64(ctr, key)threefry2x64_R(threefry2x64_rounds, ctr, key); |
428 | |
429 | rnd[0] = gaussian_table[(rand.v[0] >> 48) & GAUSS_MASK((1 << 14) - 1)]; |
430 | rnd[1] = gaussian_table[(rand.v[0] >> 32) & GAUSS_MASK((1 << 14) - 1)]; |
431 | rnd[2] = gaussian_table[(rand.v[0] >> 16) & GAUSS_MASK((1 << 14) - 1)]; |
432 | } |
433 | |
434 | void |
435 | gmx_rng_cycle_6gaussian_table(gmx_int64_t ctr1, gmx_int64_t ctr2, |
436 | gmx_int64_t key1, gmx_int64_t key2, |
437 | real* rnd) |
438 | { |
439 | threefry2x64_ctr_t ctr = {{ctr1, ctr2}}; |
440 | threefry2x64_key_t key = {{key1, key2}}; |
441 | threefry2x64_ctr_t rand = threefry2x64(ctr, key)threefry2x64_R(threefry2x64_rounds, ctr, key); |
442 | |
443 | rnd[0] = gaussian_table[(rand.v[0] >> 48) & GAUSS_MASK((1 << 14) - 1)]; |
444 | rnd[1] = gaussian_table[(rand.v[0] >> 32) & GAUSS_MASK((1 << 14) - 1)]; |
445 | rnd[2] = gaussian_table[(rand.v[0] >> 16) & GAUSS_MASK((1 << 14) - 1)]; |
446 | rnd[3] = gaussian_table[(rand.v[1] >> 48) & GAUSS_MASK((1 << 14) - 1)]; |
447 | rnd[4] = gaussian_table[(rand.v[1] >> 32) & GAUSS_MASK((1 << 14) - 1)]; |
448 | rnd[5] = gaussian_table[(rand.v[1] >> 16) & GAUSS_MASK((1 << 14) - 1)]; |
449 | } |