Manually sort some includes in src/gromacs

[alexxy/gromacs.git] / src / gromacs / random / random.c

diff --git a/src/gromacs/random/random.c b/src/gromacs/random/random.c
index 19007cf4d52b69d323d2f64a047c2e8a78e7025a..9150e64f3eb12f7425d02e5f09c05e7cf6f58838 100644 (file)
--- a/src/gromacs/random/random.c
+++ b/src/gromacs/random/random.c
@@ -34,25 +34,28 @@
  * To help us fund GROMACS development, we humbly ask that you cite
  * the research papers on the package. Check out http://www.gromacs.org.
  */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "gmxpre.h"
 
 #include "random.h"
 
+#include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
+#include <time.h>
+
+#include "config.h"
+
 #ifdef HAVE_UNISTD_H
 #include <unistd.h>
 #endif
-#include <time.h>
-#include <math.h>
 #ifdef GMX_NATIVE_WINDOWS
 #include <process.h>
 #endif
 
 #include "gromacs/math/utilities.h"
-#include "random_gausstable.h"
+#include "gromacs/random/random_gausstable.h"
+
+#include "external/Random123-1.08/include/Random123/threefry.h"
 
 #define RNG_N 624
 #define RNG_M 397
@@ -60,15 +63,15 @@
 #define RNG_UPPER_MASK 0x80000000UL /* most significant w-r bits */
 #define RNG_LOWER_MASK 0x7fffffffUL /* least significant r bits */
 
-/* Note that if you change the size of the Gaussian table you will also
- * have to generate new initialization data for the table in
+/* Note that if you change the size of the Gaussian table you will
+ * also have to generate new initialization data for the table in
  * gmx_random_gausstable.h
  *
  * A routine print_gaussian_table() is in contrib/random.c
  * for convenience - use it if you need a different size of the table.
  */
 #define GAUSS_TABLE 14 /* the size of the gauss table is 2^GAUSS_TABLE */
-#define GAUSS_SHIFT (32 - GAUSS_TABLE)
+#define GAUSS_MASK  ((1 << GAUSS_TABLE) - 1)
 
 
 struct gmx_rng {
@@ -209,7 +212,20 @@ gmx_rng_make_seed(void)
     long         my_pid;
 
 #ifdef HAVE_UNISTD_H
-    fp = fopen("/dev/random", "rb"); /* will return NULL if it is not present */
+    /* We never want Gromacs execution to halt 10-20 seconds while
+     * waiting for enough entropy in the random number generator.
+     * For this reason we should NOT use /dev/random, which will
+     * block in cases like that. That will cause all sorts of
+     * Gromacs programs to block ~20 seconds while waiting for a
+     * super-random-number to generate cool quotes. Apart from the
+     * minor irritation, it is really bad behavior of a program
+     * to abuse the system random numbers like that - other programs
+     * need them too.
+     * For this reason, we ONLY try to get random numbers from
+     * the pseudo-random stream /dev/urandom, and use other means
+     * if it is not present (in this case fopen() returns NULL).
+     */
+    fp = fopen("/dev/urandom", "rb");
 #else
     fp = NULL;
 #endif
@@ -374,10 +390,8 @@ gmx_rng_uniform_real(gmx_rng_t rng)
      * we are limited to an accuracy of 1e-7.
      */
 }
-
-
-
 real
+
 gmx_rng_gaussian_table(gmx_rng_t rng)
 {
     unsigned int i;
@@ -385,5 +399,52 @@ gmx_rng_gaussian_table(gmx_rng_t rng)
     i = gmx_rng_uniform_uint32(rng);
 
     /* The Gaussian table is a static constant in this file */
-    return gaussian_table[i >> GAUSS_SHIFT];
+    return gaussian_table[i >> (32 - GAUSS_TABLE)];
+}
+
+void
+gmx_rng_cycle_2uniform(gmx_int64_t ctr1, gmx_int64_t ctr2,
+                       gmx_int64_t key1, gmx_int64_t key2,
+                       double* rnd)
+{
+    const gmx_int64_t  mask_53bits     = 0x1FFFFFFFFFFFFF;
+    const double       two_power_min53 = 1.0/9007199254740992.0;
+
+    threefry2x64_ctr_t ctr  = {{ctr1, ctr2}};
+    threefry2x64_key_t key  = {{key1, key2}};
+    threefry2x64_ctr_t rand = threefry2x64(ctr, key);
+
+    rnd[0] = (rand.v[0] & mask_53bits)*two_power_min53;
+    rnd[1] = (rand.v[1] & mask_53bits)*two_power_min53;
+}
+
+void
+gmx_rng_cycle_3gaussian_table(gmx_int64_t ctr1, gmx_int64_t ctr2,
+                              gmx_int64_t key1, gmx_int64_t key2,
+                              real* rnd)
+{
+    threefry2x64_ctr_t ctr  = {{ctr1, ctr2}};
+    threefry2x64_key_t key  = {{key1, key2}};
+    threefry2x64_ctr_t rand = threefry2x64(ctr, key);
+
+    rnd[0] = gaussian_table[(rand.v[0] >> 48) & GAUSS_MASK];
+    rnd[1] = gaussian_table[(rand.v[0] >> 32) & GAUSS_MASK];
+    rnd[2] = gaussian_table[(rand.v[0] >> 16) & GAUSS_MASK];
+}
+
+void
+gmx_rng_cycle_6gaussian_table(gmx_int64_t ctr1, gmx_int64_t ctr2,
+                              gmx_int64_t key1, gmx_int64_t key2,
+                              real* rnd)
+{
+    threefry2x64_ctr_t ctr  = {{ctr1, ctr2}};
+    threefry2x64_key_t key  = {{key1, key2}};
+    threefry2x64_ctr_t rand = threefry2x64(ctr, key);
+
+    rnd[0] = gaussian_table[(rand.v[0] >> 48) & GAUSS_MASK];
+    rnd[1] = gaussian_table[(rand.v[0] >> 32) & GAUSS_MASK];
+    rnd[2] = gaussian_table[(rand.v[0] >> 16) & GAUSS_MASK];
+    rnd[3] = gaussian_table[(rand.v[1] >> 48) & GAUSS_MASK];
+    rnd[4] = gaussian_table[(rand.v[1] >> 32) & GAUSS_MASK];
+    rnd[5] = gaussian_table[(rand.v[1] >> 16) & GAUSS_MASK];
 }