2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2008, The GROMACS development team.
6 * Copyright (c) 2012,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
51 #ifdef GMX_NATIVE_WINDOWS
55 #include "gromacs/math/utilities.h"
56 #include "gromacs/random/random_gausstable.h"
58 #include "external/Random123-1.08/include/Random123/threefry.h"
62 #define RNG_MATRIX_A 0x9908b0dfUL /* constant vector a */
63 #define RNG_UPPER_MASK 0x80000000UL /* most significant w-r bits */
64 #define RNG_LOWER_MASK 0x7fffffffUL /* least significant r bits */
66 /* Note that if you change the size of the Gaussian table you will
67 * also have to generate new initialization data for the table in
68 * gmx_random_gausstable.h
70 * A routine print_gaussian_table() is in contrib/random.c
71 * for convenience - use it if you need a different size of the table.
73 #define GAUSS_TABLE 14 /* the size of the gauss table is 2^GAUSS_TABLE */
74 #define GAUSS_MASK ((1 << GAUSS_TABLE) - 1)
78 unsigned int mt[RNG_N];
94 gmx_rng_init(unsigned int seed)
98 if ((rng = (struct gmx_rng *)malloc(sizeof(struct gmx_rng))) == NULL)
103 rng->has_saved = 0; /* no saved gaussian number yet */
105 rng->mt[0] = seed & 0xffffffffUL;
106 for (rng->mti = 1; rng->mti < RNG_N; rng->mti++)
109 (1812433253UL * (rng->mt[rng->mti-1] ^
110 (rng->mt[rng->mti-1] >> 30)) + rng->mti);
111 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
112 /* In the previous versions, MSBs of the seed affect */
113 /* only MSBs of the array mt[]. */
114 /* 2002/01/09 modified by Makoto Matsumoto */
115 rng->mt[rng->mti] &= 0xffffffffUL;
116 /* for >32 bit machines */
122 gmx_rng_init_array(unsigned int seed[], int seed_length)
127 if ((rng = gmx_rng_init(19650218UL)) == NULL)
133 k = (RNG_N > seed_length ? RNG_N : seed_length);
136 rng->mt[i] = (rng->mt[i] ^ ((rng->mt[i-1] ^
137 (rng->mt[i-1] >> 30)) * 1664525UL))
138 + seed[j] + j; /* non linear */
139 rng->mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
143 rng->mt[0] = rng->mt[RNG_N-1]; i = 1;
145 if (j >= seed_length)
150 for (k = RNG_N-1; k; k--)
152 rng->mt[i] = (rng->mt[i] ^ ((rng->mt[i-1] ^
153 (rng->mt[i-1] >> 30)) *
155 - i; /* non linear */
156 rng->mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
160 rng->mt[0] = rng->mt[RNG_N-1]; i = 1;
164 rng->mt[0] = 0x80000000UL;
165 /* MSB is 1; assuring non-zero initial array */
171 gmx_rng_destroy(gmx_rng_t rng)
181 gmx_rng_get_state(gmx_rng_t rng, unsigned int *mt, int *mti)
185 for (i = 0; i < RNG_N; i++)
194 gmx_rng_set_state(gmx_rng_t rng, unsigned int *mt, int mti)
198 for (i = 0; i < RNG_N; i++)
207 gmx_rng_make_seed(void)
215 /* We never want Gromacs execution to halt 10-20 seconds while
216 * waiting for enough entropy in the random number generator.
217 * For this reason we should NOT use /dev/random, which will
218 * block in cases like that. That will cause all sorts of
219 * Gromacs programs to block ~20 seconds while waiting for a
220 * super-random-number to generate cool quotes. Apart from the
221 * minor irritation, it is really bad behavior of a program
222 * to abuse the system random numbers like that - other programs
224 * For this reason, we ONLY try to get random numbers from
225 * the pseudo-random stream /dev/urandom, and use other means
226 * if it is not present (in this case fopen() returns NULL).
228 fp = fopen("/dev/urandom", "rb");
234 ret = fread(&data, sizeof(unsigned int), 1, fp);
239 /* No random device available, use time-of-day and process id */
240 #ifdef GMX_NATIVE_WINDOWS
241 my_pid = (long)_getpid();
243 my_pid = (long)getpid();
245 data = (unsigned int)(((long)time(NULL)+my_pid) % (long)1000000);
251 /* The random number state contains RNG_N entries that are returned one by
252 * one as random numbers. When we run out of them, this routine is called to
253 * regenerate RNG_N new entries.
256 gmx_rng_update(gmx_rng_t rng)
258 unsigned int lastx, x1, x2, y, *mt;
260 const unsigned int mag01[2] = {0x0UL, RNG_MATRIX_A};
261 /* mag01[x] = x * MATRIX_A for x=0,1 */
263 /* update random numbers */
264 mt = rng->mt; /* pointer to array - avoid repeated dereferencing */
268 for (k = 0; k < RNG_N-RNG_M-3; k += 4)
271 y = (x1 & RNG_UPPER_MASK) | (x2 & RNG_LOWER_MASK);
272 mt[k] = mt[k+RNG_M] ^ (y >> 1) ^ mag01[y & 0x1UL];
274 y = (x2 & RNG_UPPER_MASK) | (x1 & RNG_LOWER_MASK);
275 mt[k+1] = mt[k+RNG_M+1] ^ (y >> 1) ^ mag01[y & 0x1UL];
277 y = (x1 & RNG_UPPER_MASK) | (x2 & RNG_LOWER_MASK);
278 mt[k+2] = mt[k+RNG_M+2] ^ (y >> 1) ^ mag01[y & 0x1UL];
280 y = (x2 & RNG_UPPER_MASK) | (x1 & RNG_LOWER_MASK);
281 mt[k+3] = mt[k+RNG_M+3] ^ (y >> 1) ^ mag01[y & 0x1UL];
284 y = (x1 & RNG_UPPER_MASK) | (x2 & RNG_LOWER_MASK);
285 mt[k] = mt[k+RNG_M] ^ (y >> 1) ^ mag01[y & 0x1UL];
288 y = (x2 & RNG_UPPER_MASK) | (x1 & RNG_LOWER_MASK);
289 mt[k] = mt[k+RNG_M] ^ (y >> 1) ^ mag01[y & 0x1UL];
292 y = (x1 & RNG_UPPER_MASK) | (x2 & RNG_LOWER_MASK);
293 mt[k] = mt[k+RNG_M] ^ (y >> 1) ^ mag01[y & 0x1UL];
295 for (; k < RNG_N-1; k += 4)
298 y = (x2 & RNG_UPPER_MASK) | (x1 & RNG_LOWER_MASK);
299 mt[k] = mt[k+(RNG_M-RNG_N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
301 y = (x1 & RNG_UPPER_MASK) | (x2 & RNG_LOWER_MASK);
302 mt[k+1] = mt[k+(RNG_M-RNG_N)+1] ^ (y >> 1) ^ mag01[y & 0x1UL];
304 y = (x2 & RNG_UPPER_MASK) | (x1 & RNG_LOWER_MASK);
305 mt[k+2] = mt[k+(RNG_M-RNG_N)+2] ^ (y >> 1) ^ mag01[y & 0x1UL];
307 y = (x1 & RNG_UPPER_MASK) | (x2 & RNG_LOWER_MASK);
308 mt[k+3] = mt[k+(RNG_M-RNG_N)+3] ^ (y >> 1) ^ mag01[y & 0x1UL];
310 y = (x2 & RNG_UPPER_MASK) | (mt[0] & RNG_LOWER_MASK);
311 mt[RNG_N-1] = mt[RNG_M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
318 gmx_rng_gaussian_real(gmx_rng_t rng)
325 return rng->gauss_saved;
331 x = 2.0*gmx_rng_uniform_real(rng)-1.0;
332 y = 2.0*gmx_rng_uniform_real(rng)-1.0;
335 while (r > 1.0 || r == 0.0);
337 r = sqrt(-2.0*log(r)/r);
338 rng->gauss_saved = y*r; /* save second random number */
340 return x*r; /* return first random number */
347 /* Return a random unsigned integer, i.e. 0..4294967295
348 * Provided in header file for performace reasons.
349 * Unfortunately this function cannot be inlined, since
350 * it needs to refer the internal-linkage gmx_rng_update().
353 gmx_rng_uniform_uint32(gmx_rng_t rng)
357 if (rng->mti == RNG_N)
361 y = rng->mt[rng->mti++];
364 y ^= (y << 7) & 0x9d2c5680UL;
365 y ^= (y << 15) & 0xefc60000UL;
375 /* Return a uniform floating point number on the interval 0<=x<1 */
377 gmx_rng_uniform_real(gmx_rng_t rng)
379 if (sizeof(real) == sizeof(double))
381 return ((double)gmx_rng_uniform_uint32(rng))*(1.0/4294967296.0);
385 return ((float)gmx_rng_uniform_uint32(rng))*(1.0/4294967423.0);
387 /* divided by the smallest number that will generate a
388 * single precision real number on 0<=x<1.
389 * This needs to be slightly larger than MAX_UNIT since
390 * we are limited to an accuracy of 1e-7.
395 gmx_rng_gaussian_table(gmx_rng_t rng)
399 i = gmx_rng_uniform_uint32(rng);
401 /* The Gaussian table is a static constant in this file */
402 return gaussian_table[i >> (32 - GAUSS_TABLE)];
406 gmx_rng_cycle_2uniform(gmx_int64_t ctr1, gmx_int64_t ctr2,
407 gmx_int64_t key1, gmx_int64_t key2,
410 const gmx_int64_t mask_53bits = 0x1FFFFFFFFFFFFF;
411 const double two_power_min53 = 1.0/9007199254740992.0;
413 threefry2x64_ctr_t ctr = {{ctr1, ctr2}};
414 threefry2x64_key_t key = {{key1, key2}};
415 threefry2x64_ctr_t rand = threefry2x64(ctr, key);
417 rnd[0] = (rand.v[0] & mask_53bits)*two_power_min53;
418 rnd[1] = (rand.v[1] & mask_53bits)*two_power_min53;
422 gmx_rng_cycle_3gaussian_table(gmx_int64_t ctr1, gmx_int64_t ctr2,
423 gmx_int64_t key1, gmx_int64_t key2,
426 threefry2x64_ctr_t ctr = {{ctr1, ctr2}};
427 threefry2x64_key_t key = {{key1, key2}};
428 threefry2x64_ctr_t rand = threefry2x64(ctr, key);
430 rnd[0] = gaussian_table[(rand.v[0] >> 48) & GAUSS_MASK];
431 rnd[1] = gaussian_table[(rand.v[0] >> 32) & GAUSS_MASK];
432 rnd[2] = gaussian_table[(rand.v[0] >> 16) & GAUSS_MASK];
436 gmx_rng_cycle_6gaussian_table(gmx_int64_t ctr1, gmx_int64_t ctr2,
437 gmx_int64_t key1, gmx_int64_t key2,
440 threefry2x64_ctr_t ctr = {{ctr1, ctr2}};
441 threefry2x64_key_t key = {{key1, key2}};
442 threefry2x64_ctr_t rand = threefry2x64(ctr, key);
444 rnd[0] = gaussian_table[(rand.v[0] >> 48) & GAUSS_MASK];
445 rnd[1] = gaussian_table[(rand.v[0] >> 32) & GAUSS_MASK];
446 rnd[2] = gaussian_table[(rand.v[0] >> 16) & GAUSS_MASK];
447 rnd[3] = gaussian_table[(rand.v[1] >> 48) & GAUSS_MASK];
448 rnd[4] = gaussian_table[(rand.v[1] >> 32) & GAUSS_MASK];
449 rnd[5] = gaussian_table[(rand.v[1] >> 16) & GAUSS_MASK];