5 files changed, 302 insertions, 0 deletions
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..4531139
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,17 @@
+PYTHON=python2
+
+.PHONY: all clean
+
+all: R python
+
+python:
+	$(PYTHON) setup.py build_ext --inplace
+
+R:
+	R CMD SHLIB -o tausR.so taus.c
+
+clean:
+	rm -rf build
+	rm -f taus.so
+	rm -f tausR.so
+	rm -f taus.o
diff --git a/Tausworthe.R b/Tausworthe.R
new file mode 100644
index 0000000..d35622a
--- /dev/null
+++ b/Tausworthe.R
@@ -0,0 +1,55 @@
+library(methods)
+
+dyn.load('tausR.so')
+
+TauswortheRNG <- setRefClass('TauswortheRNG',
+	fields=list(
+		    seed='numeric',
+		    state='numeric'
+		    ),
+	methods=list(
+		initialize=function(..., seed=0) {
+			seed <<- seed
+			tmp <- .Call('R_tausworthe_seed',
+			     as.integer(seed))
+			state <<- tmp[1:4]
+			callSuper(...)
+		},
+		randi=function() {
+			tmp <- .Call('R_tausworthe_rand',
+				     as.integer(state))
+			state <<- tmp[1:4]
+			return(tmp[5])
+		},
+		rand=function() {
+			r <- randi()
+			return (r * 2.3283064365387e-10)
+		}
+		)
+	)
+
+taus.seed <- function(seed=0)
+{
+	t <- TauswortheRNG(seed=seed)
+	return(t)
+}
+
+taus.rand <- function(t)
+{
+	return(t$rand())
+}
+
+taus.randi <- function(t)
+{
+	return(t$randi())
+}
+
+test.randi <- function()
+{
+	t <- TauswortheRNG(seed=112339)
+	for (i in 1:1000000) {
+		cat(t$randi(), "\n")
+	}
+}
+
+test.randi()
diff --git a/Tausworthe.py b/Tausworthe.py
new file mode 100644
index 0000000..f61fc70
--- /dev/null
+++ b/Tausworthe.py
@@ -0,0 +1,22 @@
+
+import taus
+
+class TauswortheRNG(object):
+
+    def __init__(self, seed=0):
+        self.seed = seed
+        self.state = taus.seed(seed)
+
+    def randi(self):
+        tmp = taus.rand(self.state)
+        self.state = tmp[:-1]
+        return(tmp[-1])
+
+    def rand(self):
+        return self.randi() * 2.3283064365387e-10
+
+
+if __name__ == '__main__':
+    t = TauswortheRNG(112339)
+    for i in range(1000000):
+        print(t.randi())
diff --git a/setup.py b/setup.py
new file mode 100644
index 0000000..f3c23b0
--- /dev/null
+++ b/setup.py
@@ -0,0 +1,19 @@
+
+from distutils.core import setup, Extension
+
+"""
+module1 = Extension('taus',
+        define_macros = [('TARGETPYTHON', '1')],
+        sources=['taus.c'])
+
+setup (name = 'Tausworthe RNG',
+        version = '0.1',
+        description='The Tausworthe RNG for Python and R',
+        ext_modules = [module1])
+"""
+
+setup(
+        ext_modules=[Extension("taus",
+            define_macros=[('TARGETPYTHON', '1')],
+            sources=["taus.c"])],
+        )
diff --git a/taus.c b/taus.c
new file mode 100644
index 0000000..9eba229
--- /dev/null
+++ b/taus.c
@@ -0,0 +1,189 @@
+#ifdef TARGETPYTHON
+#include "Python.h"
+#endif
+
+#ifndef TARGETPYTHON
+#define STRICT_R_HEADERS
+#include <stdio.h>
+#include <stdlib.h>
+#include <R.h>
+#include <Rinternals.h>
+#endif
+
+int lfsr113(int **state)
+{
+	unsigned long z1, z2, z3, z4, 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 & 0x7FFFFFFF;
+	(*state)[1] = z2 & 0x7FFFFFFF;
+	(*state)[2] = z3 & 0x7FFFFFFF;
+	(*state)[3] = z4 & 0x7FFFFFFF;
+	b = b & 0x7FFFFFFF;
+	
+	return(b);
+}
+
+void lfsr113_seed(unsigned long seed, int **state)
+{
+	unsigned long z1 = 2,
+		      z2 = 8,
+		      z3 = 16,
+		      z4 = 128;
+
+	z1 = (z1 * (seed + 1)) & 0x7FFFFFFF;
+	z2 = (z2 * (seed + 1)) & 0x7FFFFFFF;
+	z3 = (z3 * (seed + 1)) & 0x7FFFFFFF;
+	z4 = (z4 * (seed + 1)) & 0x7FFFFFFF;
+
+	if (*state == NULL) {
+		(*state) = malloc(sizeof(int)*4);
+	}
+
+	(*state)[0] = (int) z1;
+	(*state)[1] = (int) z2;
+	(*state)[2] = (int) z3;
+	(*state)[3] = (int) z4;
+}
+
+#ifdef TARGETPYTHON
+/*
+ *
+ * Start of Python code
+ *
+ */
+
+static char module_docstring[] = 
+"This module provides the Tausworthe RNG for R and Python simultaneously";
+
+static PyObject *taus_seed(PyObject *self, PyObject *args)
+{
+	int seed, *state = NULL;
+
+	if (!PyArg_ParseTuple(args, "i", &seed))
+		return NULL;
+
+	lfsr113_seed(seed, &state);
+	PyObject *pystate = Py_BuildValue("[i, i, i, i]",
+		       	state[0], state[1], state[2], state[3]);
+	free(state);
+	return pystate;
+}
+
+static PyObject *taus_rand(PyObject *self, PyObject *args)
+{
+	int i, value, numints, *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(int)*5);
+	numints = PyList_Size(listObj);
+	for (i=0; i<numints; i++) {
+		intObj = PyList_GetItem(listObj, i);
+		value = (int) PyLong_AsLong(intObj);
+		localstate[i] = value;
+	}
+
+	int rand = lfsr113(&localstate);
+	localstate[4] = rand;
+
+	PyObject *pystate = Py_BuildValue("[i, i, i, i, i]",
+		       	localstate[0], localstate[1], localstate[2],
+			localstate[3], rand);
+	free(localstate);
+	return pystate;
+}
+
+static PyMethodDef TausMethods[] = {
+	{"seed", taus_seed, METH_VARARGS,
+		"Seed the Tausworthe RNG."},
+	{"rand", taus_rand, METH_VARARGS,
+		"Generate a single random integer."},
+	{NULL, NULL, 0, NULL}
+};
+
+PyMODINIT_FUNC inittaus(void)
+{
+	PyObject *m = Py_InitModule3("taus", TausMethods, module_docstring);
+	if (m == NULL)
+		return;
+}
+#endif
+
+#ifndef TARGETPYTHON
+/*
+ *
+ * Start of R code
+ *
+ */
+SEXP R_tausworthe_seed(SEXP seed)
+{
+	int i, *pstate = NULL, *state = NULL;
+	int *pseed = INTEGER(seed);
+
+	SEXP Rstate = PROTECT(allocVector(INTSXP, 5));
+	pstate = INTEGER(Rstate);
+
+	lfsr113_seed(*pseed, &state);
+
+	for (i=0; i<4; i++) {
+		pstate[i] = state[i];
+	}
+	pstate[4] = -1;
+	free(state);
+
+	UNPROTECT(1);
+	return Rstate;
+}
+
+SEXP R_tausworthe_rand(SEXP state)
+{
+	int *localstate = malloc(sizeof(int)*4);
+	int *pstate = INTEGER(state);
+	int i;
+	for (i=0; i<4; i++)
+		localstate[i] = pstate[i];
+
+	int rand = lfsr113(&localstate);
+	
+	SEXP Rstate = PROTECT(allocVector(INTSXP, 5));
+	pstate = INTEGER(Rstate);
+
+	for (i=0; i<4; i++)
+		pstate[i] = localstate[i];
+	pstate[4] = rand;
+	UNPROTECT(1);
+
+	free(localstate);
+
+	return Rstate;
+}
+/*
+ *
+ * End of R code
+ *
+ */
+#endif