src/gromacs/random/random.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
  10  *
  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.
  15  *
  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.
  20  *
  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.
  25  *
  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.
  33  *
  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.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "random.h"
  40
  41 #include "config.h"
  42
  43 #include <math.h>
  44 #include <stdio.h>
  45 #include <stdlib.h>
  46 #include <time.h>
  47
  48 #include "external/Random123-1.08/include/Random123/threefry.h"
  49
  50 #include "gromacs/math/utilities.h"
  51 #include "gromacs/random/random_gausstable.h"
  52 #include "gromacs/utility/sysinfo.h"
  53
  54 #define RNG_N 624
  55 #define RNG_M 397
  56 #define RNG_MATRIX_A 0x9908b0dfUL   /* constant vector a */
  57 #define RNG_UPPER_MASK 0x80000000UL /* most significant w-r bits */
  58 #define RNG_LOWER_MASK 0x7fffffffUL /* least significant r bits */
  59
  60 /* Note that if you change the size of the Gaussian table you will
  61  * also have to generate new initialization data for the table in
  62  * gmx_random_gausstable.h
  63  *
  64  * A routine print_gaussian_table() is in contrib/random.c
  65  * for convenience - use it if you need a different size of the table.
  66  */
  67 #define GAUSS_TABLE 14 /* the size of the gauss table is 2^GAUSS_TABLE */
  68 #define GAUSS_MASK  ((1 << GAUSS_TABLE) - 1)
  69
  70
  71 struct gmx_rng {
  72     unsigned int  mt[RNG_N];
  73     int           mti;
  74     int           has_saved;
  75     double        gauss_saved;
  76 };
  77
  78
  79
  80 int
  81 gmx_rng_n(void)
  82 {
  83     return RNG_N;
  84 }
  85
  86
  87 gmx_rng_t
  88 gmx_rng_init(unsigned int seed)
  89 {
  90     struct gmx_rng *rng;
  91
  92     if ((rng = (struct gmx_rng *)malloc(sizeof(struct gmx_rng))) == NULL)
  93     {
  94         return NULL;
  95     }
  96
  97     rng->has_saved = 0; /* no saved gaussian number yet */
  98
  99     rng->mt[0] = seed & 0xffffffffUL;
 100     for (rng->mti = 1; rng->mti < RNG_N; rng->mti++)
 101     {
 102         rng->mt[rng->mti] =
 103             (1812433253UL * (rng->mt[rng->mti-1] ^
 104                              (rng->mt[rng->mti-1] >> 30)) + rng->mti);
 105         /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
 106         /* In the previous versions, MSBs of the seed affect   */
 107         /* only MSBs of the array mt[].                        */
 108         /* 2002/01/09 modified by Makoto Matsumoto             */
 109         rng->mt[rng->mti] &= 0xffffffffUL;
 110         /* for >32 bit machines */
 111     }
 112     return rng;
 113 }
 114
 115 gmx_rng_t
 116 gmx_rng_init_array(unsigned int seed[], int seed_length)
 117 {
 118     int       i, j, k;
 119     gmx_rng_t rng;
 120
 121     if ((rng = gmx_rng_init(19650218UL)) == NULL)
 122     {
 123         return NULL;
 124     }
 125
 126     i = 1; j = 0;
 127     k = (RNG_N > seed_length ? RNG_N : seed_length);
 128     for (; k; k--)
 129     {
 130         rng->mt[i] = (rng->mt[i] ^ ((rng->mt[i-1] ^
 131                                      (rng->mt[i-1] >> 30)) * 1664525UL))
 132             + seed[j] + j;          /* non linear */
 133         rng->mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
 134         i++; j++;
 135         if (i >= RNG_N)
 136         {
 137             rng->mt[0] = rng->mt[RNG_N-1]; i = 1;
 138         }
 139         if (j >= seed_length)
 140         {
 141             j = 0;
 142         }
 143     }
 144     for (k = RNG_N-1; k; k--)
 145     {
 146         rng->mt[i] = (rng->mt[i] ^ ((rng->mt[i-1] ^
 147                                      (rng->mt[i-1] >> 30)) *
 148                                     1566083941UL))
 149             - i;                    /* non linear */
 150         rng->mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
 151         i++;
 152         if (i >= RNG_N)
 153         {
 154             rng->mt[0] = rng->mt[RNG_N-1]; i = 1;
 155         }
 156     }
 157
 158     rng->mt[0] = 0x80000000UL;
 159     /* MSB is 1; assuring non-zero initial array */
 160     return rng;
 161 }
 162
 163
 164 void
 165 gmx_rng_destroy(gmx_rng_t rng)
 166 {
 167     if (rng)
 168     {
 169         free(rng);
 170     }
 171 }
 172
 173
 174 void
 175 gmx_rng_get_state(gmx_rng_t rng, unsigned int *mt, int *mti)
 176 {
 177     int i;
 178
 179     for (i = 0; i < RNG_N; i++)
 180     {
 181         mt[i] = rng->mt[i];
 182     }
 183     *mti = rng->mti;
 184 }
 185
 186
 187 void
 188 gmx_rng_set_state(gmx_rng_t rng,  unsigned int *mt, int mti)
 189 {
 190     int i;
 191
 192     for (i = 0; i < RNG_N; i++)
 193     {
 194         rng->mt[i] = mt[i];
 195     }
 196     rng->mti = mti;
 197 }
 198
 199
 200 unsigned int
 201 gmx_rng_make_seed(void)
 202 {
 203     FILE        *fp;
 204     unsigned int data;
 205     int          ret;
 206     long         my_pid;
 207
 208 #ifdef HAVE_UNISTD_H
 209     /* We never want Gromacs execution to halt 10-20 seconds while
 210      * waiting for enough entropy in the random number generator.
 211      * For this reason we should NOT use /dev/random, which will
 212      * block in cases like that. That will cause all sorts of
 213      * Gromacs programs to block ~20 seconds while waiting for a
 214      * super-random-number to generate cool quotes. Apart from the
 215      * minor irritation, it is really bad behavior of a program
 216      * to abuse the system random numbers like that - other programs
 217      * need them too.
 218      * For this reason, we ONLY try to get random numbers from
 219      * the pseudo-random stream /dev/urandom, and use other means
 220      * if it is not present (in this case fopen() returns NULL).
 221      */
 222     fp = fopen("/dev/urandom", "rb");
 223 #else
 224     fp = NULL;
 225 #endif
 226     if (fp != NULL)
 227     {
 228         ret = fread(&data, sizeof(unsigned int), 1, fp);
 229         fclose(fp);
 230     }
 231     else
 232     {
 233         /* No random device available, use time-of-day and process id */
 234         my_pid = gmx_getpid();
 235         data   = (unsigned int)(((long)time(NULL)+my_pid) % (long)1000000);
 236     }
 237     return data;
 238 }
 239
 240
 241 /* The random number state contains RNG_N entries that are returned one by
 242  * one as random numbers. When we run out of them, this routine is called to
 243  * regenerate RNG_N new entries.
 244  */
 245 static void
 246 gmx_rng_update(gmx_rng_t rng)
 247 {
 248     unsigned int       lastx, x1, x2, y, *mt;
 249     int                mti, k;
 250     const unsigned int mag01[2] = {0x0UL, RNG_MATRIX_A};
 251     /* mag01[x] = x * MATRIX_A  for x=0,1 */
 252
 253     /* update random numbers */
 254     mt  = rng->mt;  /* pointer to array - avoid repeated dereferencing */
 255     mti = rng->mti;
 256
 257     x1        = mt[0];
 258     for (k = 0; k < RNG_N-RNG_M-3; k += 4)
 259     {
 260         x2      = mt[k+1];
 261         y       = (x1 & RNG_UPPER_MASK) | (x2 & RNG_LOWER_MASK);
 262         mt[k]   = mt[k+RNG_M]   ^ (y >> 1) ^ mag01[y & 0x1UL];
 263         x1      = mt[k+2];
 264         y       = (x2 & RNG_UPPER_MASK) | (x1 & RNG_LOWER_MASK);
 265         mt[k+1] = mt[k+RNG_M+1] ^ (y >> 1) ^ mag01[y & 0x1UL];
 266         x2      = mt[k+3];
 267         y       = (x1 & RNG_UPPER_MASK) | (x2 & RNG_LOWER_MASK);
 268         mt[k+2] = mt[k+RNG_M+2] ^ (y >> 1) ^ mag01[y & 0x1UL];
 269         x1      = mt[k+4];
 270         y       = (x2 & RNG_UPPER_MASK) | (x1 & RNG_LOWER_MASK);
 271         mt[k+3] = mt[k+RNG_M+3] ^ (y >> 1) ^ mag01[y & 0x1UL];
 272     }
 273     x2        = mt[k+1];
 274     y         = (x1 & RNG_UPPER_MASK) | (x2 & RNG_LOWER_MASK);
 275     mt[k]     = mt[k+RNG_M] ^ (y >> 1) ^ mag01[y & 0x1UL];
 276     k++;
 277     x1        = mt[k+1];
 278     y         = (x2 & RNG_UPPER_MASK) | (x1 & RNG_LOWER_MASK);
 279     mt[k]     = mt[k+RNG_M] ^ (y >> 1) ^ mag01[y & 0x1UL];
 280     k++;
 281     x2        = mt[k+1];
 282     y         = (x1 & RNG_UPPER_MASK) | (x2 & RNG_LOWER_MASK);
 283     mt[k]     = mt[k+RNG_M] ^ (y >> 1) ^ mag01[y & 0x1UL];
 284     k++;
 285     for (; k < RNG_N-1; k += 4)
 286     {
 287         x1      = mt[k+1];
 288         y       = (x2 & RNG_UPPER_MASK) | (x1 & RNG_LOWER_MASK);
 289         mt[k]   = mt[k+(RNG_M-RNG_N)]   ^ (y >> 1) ^ mag01[y & 0x1UL];
 290         x2      = mt[k+2];
 291         y       = (x1 & RNG_UPPER_MASK) | (x2 & RNG_LOWER_MASK);
 292         mt[k+1] = mt[k+(RNG_M-RNG_N)+1] ^ (y >> 1) ^ mag01[y & 0x1UL];
 293         x1      = mt[k+3];
 294         y       = (x2 & RNG_UPPER_MASK) | (x1 & RNG_LOWER_MASK);
 295         mt[k+2] = mt[k+(RNG_M-RNG_N)+2] ^ (y >> 1) ^ mag01[y & 0x1UL];
 296         x2      = mt[k+4];
 297         y       = (x1 & RNG_UPPER_MASK) | (x2 & RNG_LOWER_MASK);
 298         mt[k+3] = mt[k+(RNG_M-RNG_N)+3] ^ (y >> 1) ^ mag01[y & 0x1UL];
 299     }
 300     y           = (x2 & RNG_UPPER_MASK) | (mt[0] & RNG_LOWER_MASK);
 301     mt[RNG_N-1] = mt[RNG_M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
 302
 303     rng->mti = 0;
 304 }
 305
 306
 307 real
 308 gmx_rng_gaussian_real(gmx_rng_t rng)
 309 {
 310     real x, y, r;
 311
 312     if (rng->has_saved)
 313     {
 314         rng->has_saved = 0;
 315         return rng->gauss_saved;
 316     }
 317     else
 318     {
 319         do
 320         {
 321             x = 2.0*gmx_rng_uniform_real(rng)-1.0;
 322             y = 2.0*gmx_rng_uniform_real(rng)-1.0;
 323             r = x*x+y*y;
 324         }
 325         while (r > 1.0 || r == 0.0);
 326
 327         r                = sqrt(-2.0*log(r)/r);
 328         rng->gauss_saved = y*r; /* save second random number */
 329         rng->has_saved   = 1;
 330         return x*r;             /* return first random number */
 331     }
 332 }
 333
 334
 335
 336
 337 /* Return a random unsigned integer, i.e. 0..4294967295
 338  * Provided in header file for performace reasons.
 339  * Unfortunately this function cannot be inlined, since
 340  * it needs to refer the internal-linkage gmx_rng_update().
 341  */
 342 unsigned int
 343 gmx_rng_uniform_uint32(gmx_rng_t rng)
 344 {
 345     unsigned int y;
 346
 347     if (rng->mti == RNG_N)
 348     {
 349         gmx_rng_update(rng);
 350     }
 351     y = rng->mt[rng->mti++];
 352
 353     y ^= (y >> 11);
 354     y ^= (y << 7) & 0x9d2c5680UL;
 355     y ^= (y << 15) & 0xefc60000UL;
 356     y ^= (y >> 18);
 357
 358     return y;
 359 }
 360
 361
 362
 363
 364
 365 /* Return a uniform floating point number on the interval 0<=x<1 */
 366 real
 367 gmx_rng_uniform_real(gmx_rng_t rng)
 368 {
 369     if (sizeof(real) == sizeof(double))
 370     {
 371         return ((double)gmx_rng_uniform_uint32(rng))*(1.0/4294967296.0);
 372     }
 373     else
 374     {
 375         return ((float)gmx_rng_uniform_uint32(rng))*(1.0/4294967423.0);
 376     }
 377     /* divided by the smallest number that will generate a
 378      * single precision real number on 0<=x<1.
 379      * This needs to be slightly larger than MAX_UNIT since
 380      * we are limited to an accuracy of 1e-7.
 381      */
 382 }
 383 real
 384
 385 gmx_rng_gaussian_table(gmx_rng_t rng)
 386 {
 387     unsigned int i;
 388
 389     i = gmx_rng_uniform_uint32(rng);
 390
 391     /* The Gaussian table is a static constant in this file */
 392     return gaussian_table[i >> (32 - GAUSS_TABLE)];
 393 }
 394
 395 void
 396 gmx_rng_cycle_2uniform(gmx_int64_t ctr1, gmx_int64_t ctr2,
 397                        gmx_int64_t key1, gmx_int64_t key2,
 398                        double* rnd)
 399 {
 400     const gmx_int64_t  mask_53bits     = 0x1FFFFFFFFFFFFF;
 401     const double       two_power_min53 = 1.0/9007199254740992.0;
 402
 403     threefry2x64_ctr_t ctr  = {{ctr1, ctr2}};
 404     threefry2x64_key_t key  = {{key1, key2}};
 405     threefry2x64_ctr_t rand = threefry2x64(ctr, key);
 406
 407     rnd[0] = (rand.v[0] & mask_53bits)*two_power_min53;
 408     rnd[1] = (rand.v[1] & mask_53bits)*two_power_min53;
 409 }
 410
 411 void
 412 gmx_rng_cycle_3gaussian_table(gmx_int64_t ctr1, gmx_int64_t ctr2,
 413                               gmx_int64_t key1, gmx_int64_t key2,
 414                               real* rnd)
 415 {
 416     threefry2x64_ctr_t ctr  = {{ctr1, ctr2}};
 417     threefry2x64_key_t key  = {{key1, key2}};
 418     threefry2x64_ctr_t rand = threefry2x64(ctr, key);
 419
 420     rnd[0] = gaussian_table[(rand.v[0] >> 48) & GAUSS_MASK];
 421     rnd[1] = gaussian_table[(rand.v[0] >> 32) & GAUSS_MASK];
 422     rnd[2] = gaussian_table[(rand.v[0] >> 16) & GAUSS_MASK];
 423 }
 424
 425 void
 426 gmx_rng_cycle_6gaussian_table(gmx_int64_t ctr1, gmx_int64_t ctr2,
 427                               gmx_int64_t key1, gmx_int64_t key2,
 428                               real* rnd)
 429 {
 430     threefry2x64_ctr_t ctr  = {{ctr1, ctr2}};
 431     threefry2x64_key_t key  = {{key1, key2}};
 432     threefry2x64_ctr_t rand = threefry2x64(ctr, key);
 433
 434     rnd[0] = gaussian_table[(rand.v[0] >> 48) & GAUSS_MASK];
 435     rnd[1] = gaussian_table[(rand.v[0] >> 32) & GAUSS_MASK];
 436     rnd[2] = gaussian_table[(rand.v[0] >> 16) & GAUSS_MASK];
 437     rnd[3] = gaussian_table[(rand.v[1] >> 48) & GAUSS_MASK];
 438     rnd[4] = gaussian_table[(rand.v[1] >> 32) & GAUSS_MASK];
 439     rnd[5] = gaussian_table[(rand.v[1] >> 16) & GAUSS_MASK];
 440 }