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 ab103a9d7aa4aca1ea8f0816fc0915b398ec533e..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 {
@@ -387,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;
@@ -398,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];
 }