1 files changed, 381 insertions, 0 deletions
diff --git a/new_R/src/_syncrng.c b/new_R/src/_syncrng.c
new file mode 100644
index 0000000..bd1612a
--- /dev/null
+++ b/new_R/src/_syncrng.c
@@ -0,0 +1,381 @@
+#ifdef TARGETPYTHON
+#include "Python.h"
+#include <stdint.h>
+#endif
+
+#ifndef TARGETPYTHON
+#define STRICT_R_HEADERS
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <R.h>
+#include <Rinternals.h>
+#include <R_ext/Random.h>
+#include <R_ext/Rdynload.h>
+#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;
+
+	PyObject *dblObj;
+
+	if (!PyArg_ParseTuple(args, "O", &dblObj))
+		return NULL;
+
+	seed = (uint32_t) PyLong_AsLong(dblObj);
+	lfsr113_seed(seed, &state);
+
+	PyObject *pystate = Py_BuildValue("[d, d, d, d, d]",
+			(double) state[0],
+			(double) state[1],
+			(double) state[2],
+			(double) state[3],
+			 -1.0);
+	free(state);
+	return pystate;
+}
+
+static PyObject *_syncrng_rand(PyObject *self, PyObject *args)
+{
+	int i;
+	uint32_t rand;
+	uint64_t *localstate = malloc(sizeof(uint64_t) * 4);
+
+	PyObject *listObj;
+	PyObject *dblObj;
+
+	if (!PyArg_ParseTuple(args, "O!", &PyList_Type, &listObj))
+		return NULL;
+
+	for (i=0; i<4; i++) {
+		dblObj = PyList_GetItem(listObj, i);
+		localstate[i] = (uint64_t) PyFloat_AS_DOUBLE(dblObj);
+	}
+
+	rand = lfsr113(&localstate);
+
+	PyObject *pystate = Py_BuildValue("[d, d, d, d, d]",
+			(double) localstate[0],
+			(double) localstate[1],
+			(double) localstate[2],
+			(double) localstate[3],
+			(double) rand);
+	free(localstate);
+	return pystate;
+}
+
+static PyMethodDef SyncRNGMethods[] = {
+	{"seed", _syncrng_seed, METH_VARARGS,
+		"Seed the RNG."},
+	{"rand", _syncrng_rand, METH_VARARGS,
+		"Generate a single random integer using SyncRNG."},
+	{NULL, NULL, 0, NULL}
+};
+
+#if PY_MAJOR_VERSION >= 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
+init_syncrng(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 uint32_t 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