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