aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
authorGertjan van den Burg <gertjanvandenburg@gmail.com>2021-01-14 17:51:06 +0000
committerGertjan van den Burg <gertjanvandenburg@gmail.com>2021-01-14 17:51:06 +0000
commitae2192aa2dff305ca119eddc10dfe2539852e282 (patch)
treeab8e6b4e7fd8f08aff48b19be27b50d022d49534 /R
parentPut license file back (diff)
downloadSyncRNG-ae2192aa2dff305ca119eddc10dfe2539852e282.tar.gz
SyncRNG-ae2192aa2dff305ca119eddc10dfe2539852e282.zip
Remove duplicate file from R src
Diffstat (limited to 'R')
-rw-r--r--R/src/_syncrng.c381
1 files changed, 0 insertions, 381 deletions
diff --git a/R/src/_syncrng.c b/R/src/_syncrng.c
deleted file mode 100644
index bd1612a..0000000
--- a/R/src/_syncrng.c
+++ /dev/null
@@ -1,381 +0,0 @@
-#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