Skip to content

GitLab

Commit b06fc4a7 authored by Peter Eastman's avatar Peter Eastman
Browse files

Moved SFMT into libraries folder, and restructured it to allow multiple... 

Moved SFMT into libraries folder, and restructured it to allow multiple independent random number generators
parent 625180a6
Hide whitespace changes
Inline Side-by-side
Showing with 229 additions and 150 deletions
CMakeLists.txt
View file @ b06fc4a7
...@@ -75,7 +75,7 @@ ENDIF(${CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT}) ...@@ -75,7 +75,7 @@ ENDIF(${CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT})
# The source is organized into subdirectories, but we handle them all from # The source is organized into subdirectories, but we handle them all from
# this CMakeLists file rather than letting CMake visit them as SUBDIRS. # this CMakeLists file rather than letting CMake visit them as SUBDIRS.
SET(OPENMM_SOURCE_SUBDIRS . openmmapi olla libraries/jama libraries/quern libraries/lepton platforms/reference libraries/validate) SET(OPENMM_SOURCE_SUBDIRS . openmmapi olla libraries/jama libraries/quern libraries/lepton libraries/sfmt platforms/reference libraries/validate)
# The build system will set ARCH64 for 64 bit builds, which require # The build system will set ARCH64 for 64 bit builds, which require
# use of the lib64/ library directories rather than lib/. # use of the lib64/ library directories rather than lib/.
......
platforms/reference/src/sfmt/SFMT.h → libraries/sfmt/include/sfmt/SFMT.h
View file @ b06fc4a7
/** /**
* @file SFMT.h * @file SFMT.h
* *
* @brief SIMD oriented Fast Mersenne Twister(SFMT) pseudorandom * @brief SIMD oriented Fast Mersenne Twister(SFMT) pseudorandom
* number generator * number generator
...@@ -18,7 +18,7 @@ ...@@ -18,7 +18,7 @@
* and you have to define PRIu64 and PRIx64 in this file as follows: * and you have to define PRIu64 and PRIx64 in this file as follows:
* @verbatim * @verbatim
typedef unsigned int uint32_t typedef unsigned int uint32_t
typedef unsigned long long uint64_t typedef unsigned long long uint64_t
#define PRIu64 "llu" #define PRIu64 "llu"
#define PRIx64 "llx" #define PRIx64 "llx"
@endverbatim @endverbatim
...@@ -73,89 +73,110 @@ ...@@ -73,89 +73,110 @@
#define PRE_ALWAYS inline #define PRE_ALWAYS inline
#endif #endif
uint32_t OPENMM_EXPORT gen_rand32(void); /**
uint64_t OPENMM_EXPORT gen_rand64(void); * This namespace contains the functions defined by the SFMT library.
*/
void fill_array32(uint32_t *array, int size);
void fill_array64(uint64_t *array, int size);
void OPENMM_EXPORT init_gen_rand(uint32_t seed);
void init_by_array(uint32_t *init_key, int key_length); namespace OpenMM_SFMT {
const char *get_idstring(void);
int get_min_array_size32(void); class SFMTData;
int get_min_array_size64(void); class SFMT;
OPENMM_EXPORT uint32_t gen_rand32(SFMT& sfmt);
OPENMM_EXPORT uint64_t gen_rand64(SFMT& sfmt);
OPENMM_EXPORT void fill_array32(uint32_t *array, int size, SFMT& sfmt);
OPENMM_EXPORT void fill_array64(uint64_t *array, int size, SFMT& sfmt);
OPENMM_EXPORT void init_gen_rand(uint32_t seed, SFMT& sfmt);
OPENMM_EXPORT void init_by_array(uint32_t *init_key, int key_length, SFMT& sfmt);
OPENMM_EXPORT const char* get_idstring(void);
OPENMM_EXPORT int get_min_array_size32(void);
OPENMM_EXPORT int get_min_array_size64(void);
OPENMM_EXPORT SFMTData* createSFMTData(void);
OPENMM_EXPORT void deleteSFMTData(SFMTData* data);
class SFMT {
public:
SFMT() : data(createSFMTData()) {
}
~SFMT() {
deleteSFMTData(data);
}
SFMTData* data;
};
/* These real versions are due to Isaku Wada */ /* These real versions are due to Isaku Wada */
/** generates a random number on [0,1]-real-interval */ /** generates a random number on [0,1]-real-interval */
inline static double to_real1(uint32_t v) inline static double to_real1(uint32_t v)
{ {
return v * (1.0/4294967295.0); return v * (1.0/4294967295.0);
/* divided by 2^32-1 */ /* divided by 2^32-1 */
} }
/** generates a random number on [0,1]-real-interval */ /** generates a random number on [0,1]-real-interval */
inline static double genrand_real1(void) inline static double genrand_real1(SFMT& sfmt)
{ {
return to_real1(gen_rand32()); return to_real1(gen_rand32(sfmt));
} }
/** generates a random number on [0,1)-real-interval */ /** generates a random number on [0,1)-real-interval */
inline static double to_real2(uint32_t v) inline static double to_real2(uint32_t v)
{ {
return v * (1.0/4294967296.0); return v * (1.0/4294967296.0);
/* divided by 2^32 */ /* divided by 2^32 */
} }
/** generates a random number on [0,1)-real-interval */ /** generates a random number on [0,1)-real-interval */
inline static double genrand_real2(void) inline static double genrand_real2(SFMT& sfmt)
{ {
return to_real2(gen_rand32()); return to_real2(gen_rand32(sfmt));
} }
/** generates a random number on (0,1)-real-interval */ /** generates a random number on (0,1)-real-interval */
inline static double to_real3(uint32_t v) inline static double to_real3(uint32_t v)
{ {
return (((double)v) + 0.5)*(1.0/4294967296.0); return (((double)v) + 0.5)*(1.0/4294967296.0);
/* divided by 2^32 */ /* divided by 2^32 */
} }
/** generates a random number on (0,1)-real-interval */ /** generates a random number on (0,1)-real-interval */
inline static double genrand_real3(void) inline static double genrand_real3(SFMT& sfmt)
{ {
return to_real3(gen_rand32()); return to_real3(gen_rand32(sfmt));
} }
/** These real versions are due to Isaku Wada */ /** These real versions are due to Isaku Wada */
/** generates a random number on [0,1) with 53-bit resolution*/ /** generates a random number on [0,1) with 53-bit resolution*/
inline static double to_res53(uint64_t v) inline static double to_res53(uint64_t v)
{ {
return v * (1.0/18446744073709551616.0L); return v * (1.0/18446744073709551616.0L);
} }
/** generates a random number on [0,1) with 53-bit resolution from two /** generates a random number on [0,1) with 53-bit resolution from two
* 32 bit integers */ * 32 bit integers */
inline static double to_res53_mix(uint32_t x, uint32_t y) inline static double to_res53_mix(uint32_t x, uint32_t y)
{ {
return to_res53(x | ((uint64_t)y << 32)); return to_res53(x | ((uint64_t)y << 32));
} }
/** generates a random number on [0,1) with 53-bit resolution /** generates a random number on [0,1) with 53-bit resolution
*/ */
inline static double genrand_res53(void) inline static double genrand_res53(SFMT& sfmt)
{ {
return to_res53(gen_rand64()); return to_res53(gen_rand64(sfmt));
} }
/** generates a random number on [0,1) with 53-bit resolution /** generates a random number on [0,1) with 53-bit resolution
using 32bit integer. using 32bit integer.
*/ */
inline static double genrand_res53_mix(void) inline static double genrand_res53_mix(SFMT& sfmt)
{ {
uint32_t x, y; uint32_t x, y;
x = gen_rand32(); x = gen_rand32(sfmt);
y = gen_rand32(); y = gen_rand32(sfmt);
return to_res53_mix(x, y); return to_res53_mix(x, y);
} }
} // OpenMM_SFMT
#endif #endif
platforms/reference/src/sfmt/SFMT.cpp → libraries/sfmt/src/SFMT.cpp
View file @ b06fc4a7
/** /**
* @file SFMT.c * @file SFMT.c
* @brief SIMD oriented Fast Mersenne Twister(SFMT) * @brief SIMD oriented Fast Mersenne Twister(SFMT)
* *
...@@ -10,10 +10,11 @@ ...@@ -10,10 +10,11 @@
* *
* The new BSD License is applied to this software, see LICENSE.txt * The new BSD License is applied to this software, see LICENSE.txt
*/ */
#include <string.h> #include "sfmt/SFMT.h"
#include <assert.h> #include "sfmt/SFMT-params.h"
#include "SFMT.h"
#include "SFMT-params.h" #include <cstring>
#include <cassert>
#if defined(__BIG_ENDIAN__) && !defined(__amd64) && !defined(BIG_ENDIAN64) #if defined(__BIG_ENDIAN__) && !defined(__amd64) && !defined(BIG_ENDIAN64)
#define BIG_ENDIAN64 1 #define BIG_ENDIAN64 1
...@@ -66,23 +67,55 @@ typedef struct W128_T w128_t; ...@@ -66,23 +67,55 @@ typedef struct W128_T w128_t;
/*-------------------------------------- /*--------------------------------------
FILE GLOBAL VARIABLES FILE GLOBAL VARIABLES
internal state, index counter and flag internal state, index counter and flag
--------------------------------------*/ --------------------------------------*/
/** the 128-bit internal state array */ ///** the 128-bit internal state array */
static w128_t sfmt[N]; //static w128_t sfmt[N];
/** the 32bit integer pointer to the 128-bit internal state array */ ///** the 32bit integer pointer to the 128-bit internal state array */
static uint32_t *psfmt32 = &sfmt[0].u[0]; //static uint32_t *psfmt32 = &sfmt[0].u[0];
//#if !defined(BIG_ENDIAN64) || defined(ONLY64)
///** the 64bit integer pointer to the 128-bit internal state array */
//static uint64_t *psfmt64 = (uint64_t *)&sfmt[0].u[0];
//#endif
///** index counter to the 32-bit internal state array */
//static int idx;
///** a flag: it is 0 if and only if the internal state is not yet
// * initialized. */
//static int initialized = 0;
///** a parity check vector which certificate the period of 2^{MEXP} */
//static uint32_t parity[4] = {PARITY1, PARITY2, PARITY3, PARITY4};
namespace OpenMM_SFMT {
class SFMTData {
public:
/** the 128-bit internal state array */
w128_t sfmt[N];
/** the 32bit integer pointer to the 128-bit internal state array */
uint32_t *psfmt32;
#if !defined(BIG_ENDIAN64) || defined(ONLY64)
/** the 64bit integer pointer to the 128-bit internal state array */
uint64_t *psfmt64;
#endif
/** index counter to the 32-bit internal state array */
int idx;
/** a flag: it is 0 if and only if the internal state is not yet
* initialized. */
int initialized;
/** a parity check vector which certificate the period of 2^{MEXP} */
uint32_t parity[4];
SFMTData() {
psfmt32 = &sfmt[0].u[0];
#if !defined(BIG_ENDIAN64) || defined(ONLY64) #if !defined(BIG_ENDIAN64) || defined(ONLY64)
/** the 64bit integer pointer to the 128-bit internal state array */ psfmt64 = (uint64_t *)&sfmt[0].u[0];
static uint64_t *psfmt64 = (uint64_t *)&sfmt[0].u[0];
#endif #endif
/** index counter to the 32-bit internal state array */ initialized = 0;
static int idx; parity[0] = PARITY1;
/** a flag: it is 0 if and only if the internal state is not yet parity[1] = PARITY2;
* initialized. */ parity[2] = PARITY3;
static int initialized = 0; parity[3] = PARITY4;
/** a parity check vector which certificate the period of 2^{MEXP} */ }
static uint32_t parity[4] = {PARITY1, PARITY2, PARITY3, PARITY4}; };
/*---------------- /*----------------
STATIC FUNCTIONS STATIC FUNCTIONS
...@@ -90,11 +123,11 @@ static uint32_t parity[4] = {PARITY1, PARITY2, PARITY3, PARITY4}; ...@@ -90,11 +123,11 @@ static uint32_t parity[4] = {PARITY1, PARITY2, PARITY3, PARITY4};
inline static int idxof(int i); inline static int idxof(int i);
inline static void rshift128(w128_t *out, w128_t const *in, int shift); inline static void rshift128(w128_t *out, w128_t const *in, int shift);
inline static void lshift128(w128_t *out, w128_t const *in, int shift); inline static void lshift128(w128_t *out, w128_t const *in, int shift);
inline static void gen_rand_all(void); inline static void gen_rand_all(SFMT& sfmt);
inline static void gen_rand_array(w128_t *array, int size); inline static void gen_rand_array(w128_t *array, int size, SFMT& sfmt);
inline static uint32_t func1(uint32_t x); inline static uint32_t func1(uint32_t x);
inline static uint32_t func2(uint32_t x); inline static uint32_t func2(uint32_t x);
static void period_certification(void); static void period_certification(SFMT& sfmt);
#if defined(BIG_ENDIAN64) && !defined(ONLY64) #if defined(BIG_ENDIAN64) && !defined(ONLY64)
inline static void swap(w128_t *array, int size); inline static void swap(w128_t *array, int size);
#endif #endif
...@@ -106,7 +139,7 @@ inline static void swap(w128_t *array, int size); ...@@ -106,7 +139,7 @@ inline static void swap(w128_t *array, int size);
#endif #endif
/** /**
* This function simulate a 64-bit index of LITTLE ENDIAN * This function simulate a 64-bit index of LITTLE ENDIAN
* in BIG ENDIAN machine. * in BIG ENDIAN machine.
*/ */
#ifdef ONLY64 #ifdef ONLY64
...@@ -205,7 +238,6 @@ inline static void lshift128(w128_t *out, w128_t const *in, int shift) { ...@@ -205,7 +238,6 @@ inline static void lshift128(w128_t *out, w128_t const *in, int shift) {
* @param c a 128-bit part of the internal state array * @param c a 128-bit part of the internal state array
* @param d a 128-bit part of the internal state array * @param d a 128-bit part of the internal state array
*/ */
#if (!defined(HAVE_ALTIVEC)) && (!defined(HAVE_SSE2))
#ifdef ONLY64 #ifdef ONLY64
inline static void do_recursion(w128_t *r, w128_t *a, w128_t *b, w128_t *c, inline static void do_recursion(w128_t *r, w128_t *a, w128_t *b, w128_t *c,
w128_t *d) { w128_t *d) {
...@@ -214,13 +246,13 @@ inline static void do_recursion(w128_t *r, w128_t *a, w128_t *b, w128_t *c, ...@@ -214,13 +246,13 @@ inline static void do_recursion(w128_t *r, w128_t *a, w128_t *b, w128_t *c,
lshift128(&x, a, SL2); lshift128(&x, a, SL2);
rshift128(&y, c, SR2); rshift128(&y, c, SR2);
r->u[0] = a->u[0] ^ x.u[0] ^ ((b->u[0] >> SR1) & MSK2) ^ y.u[0] r->u[0] = a->u[0] ^ x.u[0] ^ ((b->u[0] >> SR1) & MSK2) ^ y.u[0]
^ (d->u[0] << SL1); ^ (d->u[0] << SL1);
r->u[1] = a->u[1] ^ x.u[1] ^ ((b->u[1] >> SR1) & MSK1) ^ y.u[1] r->u[1] = a->u[1] ^ x.u[1] ^ ((b->u[1] >> SR1) & MSK1) ^ y.u[1]
^ (d->u[1] << SL1); ^ (d->u[1] << SL1);
r->u[2] = a->u[2] ^ x.u[2] ^ ((b->u[2] >> SR1) & MSK4) ^ y.u[2] r->u[2] = a->u[2] ^ x.u[2] ^ ((b->u[2] >> SR1) & MSK4) ^ y.u[2]
^ (d->u[2] << SL1); ^ (d->u[2] << SL1);
r->u[3] = a->u[3] ^ x.u[3] ^ ((b->u[3] >> SR1) & MSK3) ^ y.u[3] r->u[3] = a->u[3] ^ x.u[3] ^ ((b->u[3] >> SR1) & MSK3) ^ y.u[3]
^ (d->u[3] << SL1); ^ (d->u[3] << SL1);
} }
#else #else
...@@ -231,38 +263,38 @@ inline static void do_recursion(w128_t *r, w128_t *a, w128_t *b, w128_t *c, ...@@ -231,38 +263,38 @@ inline static void do_recursion(w128_t *r, w128_t *a, w128_t *b, w128_t *c,
lshift128(&x, a, SL2); lshift128(&x, a, SL2);
rshift128(&y, c, SR2); rshift128(&y, c, SR2);
r->u[0] = a->u[0] ^ x.u[0] ^ ((b->u[0] >> SR1) & MSK1) ^ y.u[0] r->u[0] = a->u[0] ^ x.u[0] ^ ((b->u[0] >> SR1) & MSK1) ^ y.u[0]
^ (d->u[0] << SL1); ^ (d->u[0] << SL1);
r->u[1] = a->u[1] ^ x.u[1] ^ ((b->u[1] >> SR1) & MSK2) ^ y.u[1] r->u[1] = a->u[1] ^ x.u[1] ^ ((b->u[1] >> SR1) & MSK2) ^ y.u[1]
^ (d->u[1] << SL1); ^ (d->u[1] << SL1);
r->u[2] = a->u[2] ^ x.u[2] ^ ((b->u[2] >> SR1) & MSK3) ^ y.u[2] r->u[2] = a->u[2] ^ x.u[2] ^ ((b->u[2] >> SR1) & MSK3) ^ y.u[2]
^ (d->u[2] << SL1); ^ (d->u[2] << SL1);
r->u[3] = a->u[3] ^ x.u[3] ^ ((b->u[3] >> SR1) & MSK4) ^ y.u[3] r->u[3] = a->u[3] ^ x.u[3] ^ ((b->u[3] >> SR1) & MSK4) ^ y.u[3]
^ (d->u[3] << SL1); ^ (d->u[3] << SL1);
} }
#endif #endif
#endif
#if (!defined(HAVE_ALTIVEC)) && (!defined(HAVE_SSE2)) #if (!defined(HAVE_ALTIVEC)) && (!defined(HAVE_SSE2))
/** /**
* This function fills the internal state array with pseudorandom * This function fills the internal state array with pseudorandom
* integers. * integers.
*/ */
inline static void gen_rand_all(void) { inline static void gen_rand_all(SFMT& sfmt) {
int i; int i;
w128_t *r1, *r2; w128_t *r1, *r2;
r1 = &sfmt[N - 2]; SFMTData& data = *sfmt.data;
r2 = &sfmt[N - 1]; r1 = &(data.sfmt[N - 2]);
r2 = &(data.sfmt[N - 1]);
for (i = 0; i < N - POS1; i++) { for (i = 0; i < N - POS1; i++) {
do_recursion(&sfmt[i], &sfmt[i], &sfmt[i + POS1], r1, r2); do_recursion(&(data.sfmt[i]), &(data.sfmt[i]), &(data.sfmt[i + POS1]), r1, r2);
r1 = r2; r1 = r2;
r2 = &sfmt[i]; r2 = &(data.sfmt[i]);
} }
for (; i < N; i++) { for (; i < N; i++) {
do_recursion(&sfmt[i], &sfmt[i], &sfmt[i + POS1 - N], r1, r2); do_recursion(&(data.sfmt[i]), &(data.sfmt[i]), &(data.sfmt[i + POS1 - N]), r1, r2);
r1 = r2; r1 = r2;
r2 = &sfmt[i]; r2 = &(data.sfmt[i]);
} }
} }
...@@ -270,22 +302,23 @@ inline static void gen_rand_all(void) { ...@@ -270,22 +302,23 @@ inline static void gen_rand_all(void) {
* This function fills the user-specified array with pseudorandom * This function fills the user-specified array with pseudorandom
* integers. * integers.
* *
* @param array an 128-bit array to be filled by pseudorandom numbers. * @param array an 128-bit array to be filled by pseudorandom numbers.
* @param size number of 128-bit pseudorandom numbers to be generated. * @param size number of 128-bit pseudorandom numbers to be generated.
*/ */
inline static void gen_rand_array(w128_t *array, int size) { inline static void gen_rand_array(w128_t *array, int size, SFMT& sfmt) {
int i, j; int i, j;
w128_t *r1, *r2; w128_t *r1, *r2;
r1 = &sfmt[N - 2]; SFMTData& data = *sfmt.data;
r2 = &sfmt[N - 1]; r1 = &(data.sfmt[N - 2]);
r2 = &(data.sfmt[N - 1]);
for (i = 0; i < N - POS1; i++) { for (i = 0; i < N - POS1; i++) {
do_recursion(&array[i], &sfmt[i], &sfmt[i + POS1], r1, r2); do_recursion(&array[i], &(data.sfmt[i]), &(data.sfmt[i + POS1]), r1, r2);
r1 = r2; r1 = r2;
r2 = &array[i]; r2 = &array[i];
} }
for (; i < N; i++) { for (; i < N; i++) {
do_recursion(&array[i], &sfmt[i], &array[i + POS1 - N], r1, r2); do_recursion(&array[i], &(data.sfmt[i]), &array[i + POS1 - N], r1, r2);
r1 = r2; r1 = r2;
r2 = &array[i]; r2 = &array[i];
} }
...@@ -295,13 +328,13 @@ inline static void gen_rand_array(w128_t *array, int size) { ...@@ -295,13 +328,13 @@ inline static void gen_rand_array(w128_t *array, int size) {
r2 = &array[i]; r2 = &array[i];
} }
for (j = 0; j < 2 * N - size; j++) { for (j = 0; j < 2 * N - size; j++) {
sfmt[j] = array[j + size - N]; data.sfmt[j] = array[j + size - N];
} }
for (; i < size; i++, j++) { for (; i < size; i++, j++) {
do_recursion(&array[i], &array[i - N], &array[i + POS1 - N], r1, r2); do_recursion(&array[i], &array[i - N], &array[i + POS1 - N], r1, r2);
r1 = r2; r1 = r2;
r2 = &array[i]; r2 = &array[i];
sfmt[j] = array[i]; data.sfmt[j] = array[i];
} }
} }
#endif #endif
...@@ -344,13 +377,14 @@ static uint32_t func2(uint32_t x) { ...@@ -344,13 +377,14 @@ static uint32_t func2(uint32_t x) {
/** /**
* This function certificate the period of 2^{MEXP} * This function certificate the period of 2^{MEXP}
*/ */
static void period_certification(void) { static void period_certification(SFMT& sfmt) {
int inner = 0; int inner = 0;
int i, j; int i, j;
uint32_t work; uint32_t work;
SFMTData& data = *sfmt.data;
for (i = 0; i < 4; i++) for (i = 0; i < 4; i++)
inner ^= psfmt32[idxof(i)] & parity[i]; inner ^= data.psfmt32[idxof(i)] & data.parity[i];
for (i = 16; i > 0; i >>= 1) for (i = 16; i > 0; i >>= 1)
inner ^= inner >> i; inner ^= inner >> i;
inner &= 1; inner &= 1;
...@@ -362,8 +396,8 @@ static void period_certification(void) { ...@@ -362,8 +396,8 @@ static void period_certification(void) {
for (i = 0; i < 4; i++) { for (i = 0; i < 4; i++) {
work = 1; work = 1;
for (j = 0; j < 32; j++) { for (j = 0; j < 32; j++) {
if ((work & parity[i]) != 0) { if ((work & data.parity[i]) != 0) {
psfmt32[idxof(i)] ^= work; data.psfmt32[idxof(i)] ^= work;
return; return;
} }
work = work << 1; work = work << 1;
...@@ -407,15 +441,16 @@ int get_min_array_size64(void) { ...@@ -407,15 +441,16 @@ int get_min_array_size64(void) {
* init_gen_rand or init_by_array must be called before this function. * init_gen_rand or init_by_array must be called before this function.
* @return 32-bit pseudorandom number * @return 32-bit pseudorandom number
*/ */
uint32_t gen_rand32(void) { uint32_t gen_rand32(SFMT& sfmt) {
uint32_t r; uint32_t r;
SFMTData& data = *sfmt.data;
assert(initialized); assert(data.initialized);
if (idx >= N32) { if (data.idx >= N32) {
gen_rand_all(); gen_rand_all(sfmt);
idx = 0; data.idx = 0;
} }
r = psfmt32[idx++]; r = data.psfmt32[data.idx++];
return r; return r;
} }
#endif #endif
...@@ -423,31 +458,32 @@ uint32_t gen_rand32(void) { ...@@ -423,31 +458,32 @@ uint32_t gen_rand32(void) {
* This function generates and returns 64-bit pseudorandom number. * This function generates and returns 64-bit pseudorandom number.
* init_gen_rand or init_by_array must be called before this function. * init_gen_rand or init_by_array must be called before this function.
* The function gen_rand64 should not be called after gen_rand32, * The function gen_rand64 should not be called after gen_rand32,
* unless an initialization is again executed. * unless an initialization is again executed.
* @return 64-bit pseudorandom number * @return 64-bit pseudorandom number
*/ */
uint64_t gen_rand64(void) { uint64_t gen_rand64(SFMT& sfmt) {
#if defined(BIG_ENDIAN64) && !defined(ONLY64) #if defined(BIG_ENDIAN64) && !defined(ONLY64)
uint32_t r1, r2; uint32_t r1, r2;
#else #else
uint64_t r; uint64_t r;
#endif #endif
SFMTData& data = *sfmt.data;
assert(initialized); assert(data.initialized);
assert(idx % 2 == 0); assert(data.idx % 2 == 0);
if (idx >= N32) { if (data.idx >= N32) {
gen_rand_all(); gen_rand_all(sfmt);
idx = 0; data.idx = 0;
} }
#if defined(BIG_ENDIAN64) && !defined(ONLY64) #if defined(BIG_ENDIAN64) && !defined(ONLY64)
r1 = psfmt32[idx]; r1 = data.psfmt32[data.idx];
r2 = psfmt32[idx + 1]; r2 = data.psfmt32[data.idx + 1];
idx += 2; data.idx += 2;
return ((uint64_t)r2 << 32) | r1; return ((uint64_t)r2 << 32) | r1;
#else #else
r = psfmt64[idx / 2]; r = data.psfmt64[data.idx / 2];
idx += 2; data.idx += 2;
return r; return r;
#endif #endif
} }
...@@ -478,14 +514,15 @@ uint64_t gen_rand64(void) { ...@@ -478,14 +514,15 @@ uint64_t gen_rand64(void) {
* memory. Mac OSX doesn't have these functions, but \b malloc of OSX * memory. Mac OSX doesn't have these functions, but \b malloc of OSX
* returns the pointer to the aligned memory block. * returns the pointer to the aligned memory block.
*/ */
void fill_array32(uint32_t *array, int size) { void fill_array32(uint32_t *array, int size, SFMT& sfmt) {
assert(initialized); SFMTData& data = *sfmt.data;
assert(idx == N32); assert(data.initialized);
assert(data.idx == N32);
assert(size % 4 == 0); assert(size % 4 == 0);
assert(size >= N32); assert(size >= N32);
gen_rand_array((w128_t *)array, size / 4); gen_rand_array((w128_t *)array, size / 4, sfmt);
idx = N32; data.idx = N32;
} }
#endif #endif
...@@ -514,14 +551,15 @@ void fill_array32(uint32_t *array, int size) { ...@@ -514,14 +551,15 @@ void fill_array32(uint32_t *array, int size) {
* memory. Mac OSX doesn't have these functions, but \b malloc of OSX * memory. Mac OSX doesn't have these functions, but \b malloc of OSX
* returns the pointer to the aligned memory block. * returns the pointer to the aligned memory block.
*/ */
void fill_array64(uint64_t *array, int size) { void fill_array64(uint64_t *array, int size, SFMT& sfmt) {
assert(initialized); SFMTData& data = *sfmt.data;
assert(idx == N32); assert(data.initialized);
assert(data.idx == N32);
assert(size % 2 == 0); assert(size % 2 == 0);
assert(size >= N64); assert(size >= N64);
gen_rand_array((w128_t *)array, size / 2); gen_rand_array((w128_t *)array, size / 2, sfmt);
idx = N32; data.idx = N32;
#if defined(BIG_ENDIAN64) && !defined(ONLY64) #if defined(BIG_ENDIAN64) && !defined(ONLY64)
swap((w128_t *)array, size /2); swap((w128_t *)array, size /2);
...@@ -534,18 +572,19 @@ void fill_array64(uint64_t *array, int size) { ...@@ -534,18 +572,19 @@ void fill_array64(uint64_t *array, int size) {
* *
* @param seed a 32-bit integer used as the seed. * @param seed a 32-bit integer used as the seed.
*/ */
void init_gen_rand(uint32_t seed) { void init_gen_rand(uint32_t seed, SFMT& sfmt) {
int i; int i;
SFMTData& data = *sfmt.data;
psfmt32[idxof(0)] = seed; data.psfmt32[idxof(0)] = seed;
for (i = 1; i < N32; i++) { for (i = 1; i < N32; i++) {
psfmt32[idxof(i)] = 1812433253UL * (psfmt32[idxof(i - 1)] data.psfmt32[idxof(i)] = 1812433253UL * (data.psfmt32[idxof(i - 1)]
^ (psfmt32[idxof(i - 1)] >> 30)) ^ (data.psfmt32[idxof(i - 1)] >> 30))
+ i; + i;
} }
idx = N32; data.idx = N32;
period_certification(); period_certification(sfmt);
initialized = 1; data.initialized = 1;
} }
/** /**
...@@ -554,7 +593,7 @@ void init_gen_rand(uint32_t seed) { ...@@ -554,7 +593,7 @@ void init_gen_rand(uint32_t seed) {
* @param init_key the array of 32-bit integers, used as a seed. * @param init_key the array of 32-bit integers, used as a seed.
* @param key_length the length of init_key. * @param key_length the length of init_key.
*/ */
void init_by_array(uint32_t *init_key, int key_length) { void init_by_array(uint32_t *init_key, int key_length, SFMT& sfmt) {
int i, j, count; int i, j, count;
uint32_t r; uint32_t r;
int lag; int lag;
...@@ -572,49 +611,68 @@ void init_by_array(uint32_t *init_key, int key_length) { ...@@ -572,49 +611,68 @@ void init_by_array(uint32_t *init_key, int key_length) {
} }
mid = (size - lag) / 2; mid = (size - lag) / 2;
memset(sfmt, 0x8b, sizeof(sfmt)); SFMTData& data = *sfmt.data;
memset(data.sfmt, 0x8b, sizeof(data.sfmt));
if (key_length + 1 > N32) { if (key_length + 1 > N32) {
count = key_length + 1; count = key_length + 1;
} else { } else {
count = N32; count = N32;
} }
r = func1(psfmt32[idxof(0)] ^ psfmt32[idxof(mid)] r = func1(data.psfmt32[idxof(0)] ^ data.psfmt32[idxof(mid)]
^ psfmt32[idxof(N32 - 1)]); ^ data.psfmt32[idxof(N32 - 1)]);
psfmt32[idxof(mid)] += r; data.psfmt32[idxof(mid)] += r;
r += key_length; r += key_length;
psfmt32[idxof(mid + lag)] += r; data.psfmt32[idxof(mid + lag)] += r;
psfmt32[idxof(0)] = r; data.psfmt32[idxof(0)] = r;
count--; count--;
for (i = 1, j = 0; (j < count) && (j < key_length); j++) { for (i = 1, j = 0; (j < count) && (j < key_length); j++) {
r = func1(psfmt32[idxof(i)] ^ psfmt32[idxof((i + mid) % N32)] r = func1(data.psfmt32[idxof(i)] ^ data.psfmt32[idxof((i + mid) % N32)]
^ psfmt32[idxof((i + N32 - 1) % N32)]); ^ data.psfmt32[idxof((i + N32 - 1) % N32)]);
psfmt32[idxof((i + mid) % N32)] += r; data.psfmt32[idxof((i + mid) % N32)] += r;
r += init_key[j] + i; r += init_key[j] + i;
psfmt32[idxof((i + mid + lag) % N32)] += r; data.psfmt32[idxof((i + mid + lag) % N32)] += r;
psfmt32[idxof(i)] = r; data.psfmt32[idxof(i)] = r;
i = (i + 1) % N32; i = (i + 1) % N32;
} }
for (; j < count; j++) { for (; j < count; j++) {
r = func1(psfmt32[idxof(i)] ^ psfmt32[idxof((i + mid) % N32)] r = func1(data.psfmt32[idxof(i)] ^ data.psfmt32[idxof((i + mid) % N32)]
^ psfmt32[idxof((i + N32 - 1) % N32)]); ^ data.psfmt32[idxof((i + N32 - 1) % N32)]);
psfmt32[idxof((i + mid) % N32)] += r; data.psfmt32[idxof((i + mid) % N32)] += r;
r += i; r += i;
psfmt32[idxof((i + mid + lag) % N32)] += r; data.psfmt32[idxof((i + mid + lag) % N32)] += r;
psfmt32[idxof(i)] = r; data.psfmt32[idxof(i)] = r;
i = (i + 1) % N32; i = (i + 1) % N32;
} }
for (j = 0; j < N32; j++) { for (j = 0; j < N32; j++) {
r = func2(psfmt32[idxof(i)] + psfmt32[idxof((i + mid) % N32)] r = func2(data.psfmt32[idxof(i)] + data.psfmt32[idxof((i + mid) % N32)]
+ psfmt32[idxof((i + N32 - 1) % N32)]); + data.psfmt32[idxof((i + N32 - 1) % N32)]);
psfmt32[idxof((i + mid) % N32)] ^= r; data.psfmt32[idxof((i + mid) % N32)] ^= r;
r -= i; r -= i;
psfmt32[idxof((i + mid + lag) % N32)] ^= r; data.psfmt32[idxof((i + mid + lag) % N32)] ^= r;
psfmt32[idxof(i)] = r; data.psfmt32[idxof(i)] = r;
i = (i + 1) % N32; i = (i + 1) % N32;
} }
idx = N32; data.idx = N32;
period_certification(); period_certification(sfmt);
initialized = 1; data.initialized = 1;
} }
/**
* Create an SFMTData object. This allows outside code to create these objects without needing to know their definition.
*/
SFMTData* createSFMTData(void) {
return new SFMTData();
}
/**
* Delete an SFMTData object that was created with createSFMTData().
*/
void deleteSFMTData(SFMTData* data) {
delete data;
}
} // OpenMM_SFMT
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment