From 0b05168b80f627e175edee95a79cd0fd45f01338 Mon Sep 17 00:00:00 2001 From: Gertjan van den Burg Date: Fri, 15 Dec 2017 15:46:06 -0500 Subject: Make SyncRNG an actual user-defined RNG --- src/syncrng.c | 111 ++++++++++++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 97 insertions(+), 14 deletions(-) (limited to 'src/syncrng.c') diff --git a/src/syncrng.c b/src/syncrng.c index 025b1ca..8957047 100644 --- a/src/syncrng.c +++ b/src/syncrng.c @@ -10,6 +10,7 @@ #include #include #include +#include #endif /** @@ -36,9 +37,10 @@ * * @return a generated random number */ -uint64_t lfsr113(uint64_t **state) +uint32_t lfsr113(uint64_t **state) { - uint64_t z1, z2, z3, z4, b; + uint64_t z1, z2, z3, z4; + uint64_t b; z1 = (*state)[0]; z2 = (*state)[1]; @@ -66,7 +68,7 @@ uint64_t lfsr113(uint64_t **state) b = b & 0xFFFFFFFF; - return(b); + return((uint32_t) b); } /** @@ -80,12 +82,12 @@ uint64_t lfsr113(uint64_t **state) * @param[in] seed user supplied seed value for the RNG * @param[out] state state of the RNG */ -void lfsr113_seed(uint64_t seed, uint64_t **state) +void lfsr113_seed(uint32_t seed, uint64_t **state) { uint64_t z1 = 2, - z2 = 8, - z3 = 16, - z4 = 128; + z2 = 8, + z3 = 16, + z4 = 128; z1 = (z1 * (seed + 1)); z2 = (z2 * (seed + 1)); @@ -116,7 +118,7 @@ void lfsr113_seed(uint64_t seed, uint64_t **state) static PyObject *syncrng_seed(PyObject *self, PyObject *args) { - uint64_t seed, *state = NULL; + uint32_t seed, *state = NULL; if (!PyArg_ParseTuple(args, "k", &seed)) return NULL; @@ -131,7 +133,7 @@ static PyObject *syncrng_seed(PyObject *self, PyObject *args) static PyObject *syncrng_rand(PyObject *self, PyObject *args) { - uint64_t i, value, numints, *localstate; + uint32_t i, value, numints, *localstate; PyObject *listObj; PyObject *intObj; @@ -140,15 +142,15 @@ static PyObject *syncrng_rand(PyObject *self, PyObject *args) return NULL; // we're just assuming you would never pass more than 4 values - localstate = malloc(sizeof(uint64_t)*5); + localstate = malloc(sizeof(uint32_t)*5); numints = PyList_Size(listObj); for (i=0; i 0.259/u + 0.35) + continue; + } while (z > -log(u)); + global_R_result_normal = x; + return &global_R_result_normal; +} + /* * * End of R code -- cgit v1.2.3 From db0755d027f0f88568909e16b715c8e769333390 Mon Sep 17 00:00:00 2001 From: Gertjan van den Burg Date: Fri, 15 Dec 2017 17:25:50 -0500 Subject: fix warnings from R about routine registration --- src/syncrng.c | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) (limited to 'src/syncrng.c') diff --git a/src/syncrng.c b/src/syncrng.c index 8957047..8a0b5b4 100644 --- a/src/syncrng.c +++ b/src/syncrng.c @@ -11,6 +11,7 @@ #include #include #include +#include #endif /** @@ -221,6 +222,24 @@ initsyncrng(void) * */ +SEXP R_syncrng_seed(SEXP seed); +SEXP R_syncrng_rand(SEXP state); + +R_CallMethodDef callMethods[] = { + {"R_syncrng_seed", (DL_FUNC) &R_syncrng_seed, 1}, + {"R_syncrng_rand", (DL_FUNC) &R_syncrng_seed, 1}, + {NULL, NULL, 0} +}; +R_CMethodDef cMethods[] = { + {NULL, NULL, 0} +}; + +void R_init_myLib(DllInfo *info) +{ + R_registerRoutines(info, cMethods, callMethods, NULL, NULL); + R_useDynamicSymbols(info, TRUE); +} + // Set the seed for the generator from the reference class SEXP R_syncrng_seed(SEXP seed) { -- cgit v1.2.3 From 629af12279a9f6244ecbcb3d0021df7bc94188a3 Mon Sep 17 00:00:00 2001 From: Gertjan van den Burg Date: Fri, 15 Dec 2017 17:41:51 -0500 Subject: merge fixes from python branch --- src/syncrng.c | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) (limited to 'src/syncrng.c') diff --git a/src/syncrng.c b/src/syncrng.c index 8a0b5b4..3cc77fa 100644 --- a/src/syncrng.c +++ b/src/syncrng.c @@ -119,7 +119,8 @@ void lfsr113_seed(uint32_t seed, uint64_t **state) static PyObject *syncrng_seed(PyObject *self, PyObject *args) { - uint32_t seed, *state = NULL; + uint32_t seed; + uint64_t *state = NULL; if (!PyArg_ParseTuple(args, "k", &seed)) return NULL; @@ -134,7 +135,8 @@ static PyObject *syncrng_seed(PyObject *self, PyObject *args) static PyObject *syncrng_rand(PyObject *self, PyObject *args) { - uint32_t i, value, numints, *localstate; + uint32_t i, value, numints; + uint64_t *localstate; PyObject *listObj; PyObject *intObj; -- cgit v1.2.3 From 8f4082e9be91169a04252e77206ae15ea8554b6b Mon Sep 17 00:00:00 2001 From: Gertjan van den Burg Date: Sun, 17 Dec 2017 09:19:11 -0500 Subject: change type of global_R_nseed to match function --- src/syncrng.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src/syncrng.c') diff --git a/src/syncrng.c b/src/syncrng.c index 3cc77fa..71fe58e 100644 --- a/src/syncrng.c +++ b/src/syncrng.c @@ -298,7 +298,7 @@ SEXP R_syncrng_rand(SEXP state) */ static uint32_t global_R_seed; -static uint32_t global_R_nseed = 1; +static int global_R_nseed = 1; static double global_R_result_uniform; static double global_R_result_normal; static uint64_t *global_R_state = NULL; -- cgit v1.2.3 From 4214782b1e06ce8510d70cc2fc36a3b544508948 Mon Sep 17 00:00:00 2001 From: Gertjan van den Burg Date: Thu, 14 Jan 2021 17:15:44 +0000 Subject: Move R package to separate directory Rearranging repository --- src/syncrng.c | 375 ---------------------------------------------------------- 1 file changed, 375 deletions(-) delete mode 100644 src/syncrng.c (limited to 'src/syncrng.c') diff --git a/src/syncrng.c b/src/syncrng.c deleted file mode 100644 index 71fe58e..0000000 --- a/src/syncrng.c +++ /dev/null @@ -1,375 +0,0 @@ -#ifdef TARGETPYTHON -#include "Python.h" -#include -#endif - -#ifndef TARGETPYTHON -#define STRICT_R_HEADERS -#include -#include -#include -#include -#include -#include -#include -#endif - -/** - * @brief Generate a single random number using the capped Tausworthe RNG - * - * @details - * This generates random numbers according to the process described in [1]. As - * an additional step, the resulting random number is capped to 0xFFFFFFFF - * using a bitwise and. This is done to yield the range [0, 2^32-1]. On - * return, the state variables are updated. - * - * [1]: @article{l1996maximally, - * title={Maximally equidistributed combined Tausworthe generators}, - * author={L’ecuyer, Pierre}, - * journal={Mathematics of Computation of the American Mathematical - * Society}, - * volume={65}, - * number={213}, - * pages={203--213}, - * year={1996} - * } - * - * @param[in,out] state pointer to current state array - * - * @return a generated random number - */ -uint32_t lfsr113(uint64_t **state) -{ - uint64_t z1, z2, z3, z4; - uint64_t b; - - z1 = (*state)[0]; - z2 = (*state)[1]; - z3 = (*state)[2]; - z4 = (*state)[3]; - - b = (((z1 << 6) ^ z1) >> 13); - z1 = (((z1 & 4294967294) << 18) ^ b); - - b = (((z2 << 2) ^ z2) >> 27); - z2 = (((z2 & 4294967288) << 2) ^ b); - - b = (((z3 << 13) ^ z3) >> 21); - z3 = (((z3 & 4294967280) << 7) ^ b); - - b = (((z4 << 3) ^ z4) >> 12); - z4 = (((z4 & 4294967168) << 13) ^ b); - - b = (z1 ^ z2 ^ z3 ^ z4); - - (*state)[0] = z1; - (*state)[1] = z2; - (*state)[2] = z3; - (*state)[3] = z4; - - b = b & 0xFFFFFFFF; - - return((uint32_t) b); -} - -/** - * @brief Seed the Tausworthe RNG using a seed value - * - * @details - * This function seeds the state array using a supplied seed value. As noted - * in [1] (see lfsr113()), the values of z1, z2, z3, and z4 should be larger - * than 1, 7, 15, and 127 respectively. - * - * @param[in] seed user supplied seed value for the RNG - * @param[out] state state of the RNG - */ -void lfsr113_seed(uint32_t seed, uint64_t **state) -{ - uint64_t z1 = 2, - z2 = 8, - z3 = 16, - z4 = 128; - - z1 = (z1 * (seed + 1)); - z2 = (z2 * (seed + 1)); - z3 = (z3 * (seed + 1)); - z4 = (z4 * (seed + 1)); - - z1 = (z1 > 1) ? z1 : z1 + 1; - z2 = (z2 > 7) ? z2 : z2 + 7; - z3 = (z3 > 15) ? z3 : z3 + 15; - z4 = (z4 > 127) ? z4 : z4 + 127; - - if (*state == NULL) { - (*state) = malloc(sizeof(uint64_t)*4); - } - - (*state)[0] = z1; - (*state)[1] = z2; - (*state)[2] = z3; - (*state)[3] = z4; -} - -#ifdef TARGETPYTHON -/* - * - * Start of Python code - * - */ - -static PyObject *syncrng_seed(PyObject *self, PyObject *args) -{ - uint32_t seed; - uint64_t *state = NULL; - - if (!PyArg_ParseTuple(args, "k", &seed)) - return NULL; - - lfsr113_seed(seed, &state); - - PyObject *pystate = Py_BuildValue("[k, k, k, k]", - state[0], state[1], state[2], state[3]); - free(state); - return pystate; -} - -static PyObject *syncrng_rand(PyObject *self, PyObject *args) -{ - uint32_t i, value, numints; - uint64_t *localstate; - - PyObject *listObj; - PyObject *intObj; - - if (!PyArg_ParseTuple(args, "O!", &PyList_Type, &listObj)) - return NULL; - - // we're just assuming you would never pass more than 4 values - localstate = malloc(sizeof(uint32_t)*5); - numints = PyList_Size(listObj); - for (i=0; i= 3 - static struct PyModuleDef moduledef = { - PyModuleDef_HEAD_INIT, - "syncrng", - "Python interface to SyncRNG", - -1, - SyncRNGMethods, - NULL, - NULL, - NULL, - NULL - }; -#endif - - -static PyObject * -moduleinit(void) -{ - PyObject *m; - - #if PY_MAJOR_VERSION >= 3 - m = PyModule_Create(&moduledef); - #else - m = Py_InitModule3("syncrng", SyncRNGMethods, - "Python interface to SyncRNG"); - #endif - - return m; -} - -#if PY_MAJOR_VERSION >= 3 -PyMODINIT_FUNC -PyInit_syncrng(void) -{ - return moduleinit(); -} -#else -PyMODINIT_FUNC -initsyncrng(void) -{ - moduleinit(); -} -#endif -#endif - -#ifndef TARGETPYTHON - -/* - * - * Start of R code - * - */ - -SEXP R_syncrng_seed(SEXP seed); -SEXP R_syncrng_rand(SEXP state); - -R_CallMethodDef callMethods[] = { - {"R_syncrng_seed", (DL_FUNC) &R_syncrng_seed, 1}, - {"R_syncrng_rand", (DL_FUNC) &R_syncrng_seed, 1}, - {NULL, NULL, 0} -}; -R_CMethodDef cMethods[] = { - {NULL, NULL, 0} -}; - -void R_init_myLib(DllInfo *info) -{ - R_registerRoutines(info, cMethods, callMethods, NULL, NULL); - R_useDynamicSymbols(info, TRUE); -} - -// Set the seed for the generator from the reference class -SEXP R_syncrng_seed(SEXP seed) -{ - int i; - double *pseed = REAL(seed), - *pstate = NULL; - uint32_t useed = (uint32_t) *pseed; - uint64_t *state = NULL; - - lfsr113_seed(useed, &state); - - SEXP Rstate = PROTECT(allocVector(REALSXP, 5)); - pstate = REAL(Rstate); - for (i=0; i<4; i++) { - pstate[i] = (double) state[i]; - } - pstate[4] = -1.0; - free(state); - - UNPROTECT(1); - return Rstate; -} - -// get a random number from the reference class -SEXP R_syncrng_rand(SEXP state) -{ - uint64_t *localstate = malloc(sizeof(uint64_t)*4); - double *pstate = REAL(state); - int i; - for (i=0; i<4; i++) { - localstate[i] = (uint64_t) pstate[i]; - } - - uint32_t rand = lfsr113(&localstate); - - SEXP Rstate = PROTECT(allocVector(REALSXP, 5)); - pstate = REAL(Rstate); - - for (i=0; i<4; i++) { - pstate[i] = (double) localstate[i]; - } - pstate[4] = (double) rand; - UNPROTECT(1); - - free(localstate); - - return Rstate; -} - -/* - * The following code is used to make SyncRNG a real "user-defined" RNG - * follwing .Random.user documentation. - * - */ - -static uint32_t global_R_seed; -static int global_R_nseed = 1; -static double global_R_result_uniform; -static double global_R_result_normal; -static uint64_t *global_R_state = NULL; - -double *user_unif_rand() -{ - if (global_R_state == NULL) { - // if it's not seeded yet we seed it with 0 - global_R_seed = 0; - lfsr113_seed(global_R_seed, &global_R_state); - } - - uint32_t rand = lfsr113(&global_R_state); - global_R_result_uniform = rand * 2.3283064365387e-10; - return &global_R_result_uniform; -} - -// see: https://stackoverflow.com/q/47824450/1154005 -Int32 _unscramble(Int32 scram) -{ - int j; - for (j=0; j<50; j++) { - scram = ((scram - 1) * 2783094533); - } - return scram; -} - -// note that Int32 is "unsigned int" which is not necessarily 32 bit -void user_unif_init(Int32 seed_in) -{ - global_R_seed = seed_in; - uint32_t useed = _unscramble(seed_in); - - // destroy the previous state, we're reseeding the RNG - if (global_R_state != NULL) { - free(global_R_state); - global_R_state = NULL; - } - lfsr113_seed(useed, &global_R_state); -} - -int *user_unif_nseed() -{ - return &global_R_nseed; -} - -int *user_unif_seedloc() -{ - return (int *) &global_R_seed; -} - -double *user_norm_rand() -{ - double u, v, z, x; - do { - u = *user_unif_rand(); - v = 0.857764 * (2. * (*user_unif_rand()) - 1); - x = v/u; - z = 0.25 * x * x; - if (z < 1. - u) - break; - if (z > 0.259/u + 0.35) - continue; - } while (z > -log(u)); - global_R_result_normal = x; - return &global_R_result_normal; -} - -/* - * - * End of R code - * - */ -#endif -- cgit v1.2.3