ENH: Add support for xoroshiro128

Clean up splitmix64 to use external functions Add xoroshiro128 to show use of additional PRNG
numpy · seberg · May 28, 2019 · Feb 20, 2018 · Feb 22, 2018 · Feb 23, 2018
commit fa3bef5639b2f3223eeee7c50155824b8df991e0
diff --git a/_randomgen/README.md b/_randomgen/README.md
@@ -12,10 +12,19 @@ python setup.py develop
 ```
 
 ```ipython
-In [13]: import core_prng.generator
+In [1]: import core_prng.generator
 
-In [14]: rg = core_prng.generator.RandomGenerator()
+# Default generator is Splitmix64
+In [2]: rg = core_prng.generator.RandomGenerator()
 
-In [15]: rg.random_integer()
-Out[15]: 872337561037043212
+In [3]: rg.random_integer()
+Out[3]: 872337561037043212
+
+In [4]: from core_prng.xoroshiro128 import Xoroshiro128
+
+# Swap the generator
+In [5]: rg = core_prng.generator.RandomGenerator(Xoroshiro128())
+
+In [6]: rg.random_integer()
+Out[6]: 13370384800127340062
 ```
diff --git a/_randomgen/core_prng/common.c b/_randomgen/core_prng/common.c
diff --git a/_randomgen/core_prng/common.h b/_randomgen/core_prng/common.h
diff --git a/_randomgen/core_prng/generator.pyx b/_randomgen/core_prng/generator.pyx
@@ -3,7 +3,7 @@ cimport numpy as np
 from cpython.pycapsule cimport PyCapsule_IsValid, PyCapsule_GetPointer
 from common cimport *
 
-from core_prng.core_prng import CorePRNG
+from core_prng.splitmix64 import SplitMix64
 
 np.import_array()
 
@@ -24,18 +24,24 @@ cdef class RandomGenerator:
     """
     cdef public object __core_prng
     cdef anon_func_state anon_rng_func_state
+    cdef random_uint64_anon next_uint64
+    cdef void *rng_state
 
     def __init__(self, prng=None):
         if prng is None:
-            prng = CorePRNG()
+            prng = SplitMix64()
         self.__core_prng = prng
 
         capsule = prng._anon_func_state
         cdef const char *anon_name = "Anon CorePRNG func_state"
         if not PyCapsule_IsValid(capsule, anon_name):
             raise ValueError("Invalid pointer to anon_func_state")
         self.anon_rng_func_state = (<anon_func_state *>PyCapsule_GetPointer(capsule, anon_name))[0]
+        self.next_uint64 = <random_uint64_anon>self.anon_rng_func_state.f
+        self.rng_state = self.anon_rng_func_state.state
 
     def random_integer(self):
-        cdef random_uint64_anon f = <random_uint64_anon>self.anon_rng_func_state.f
-        return f(self.anon_rng_func_state.state)
+        return self.next_uint64(self.rng_state)
+
+    def random_double(self):
+        return (self.next_uint64(self.rng_state) >> 11) * (1.0 / 9007199254740992.0)
diff --git a/_randomgen/core_prng/core_prng.pyx → _randomgen/core_prng/splitmix64.pyx b/_randomgen/core_prng/core_prng.pyx → _randomgen/core_prng/splitmix64.pyx
@@ -6,36 +6,22 @@ from common cimport *
 
 np.import_array()
 
-cdef struct state:
-    uint64_t state
+cdef extern from "src/splitmix64/splitmix64.h":
 
-ctypedef state state_t
+    cdef struct s_splitmix64_state:
+        uint64_t state
 
-ctypedef uint64_t (*random_uint64)(state_t *st)
+    ctypedef s_splitmix64_state splitmix64_state
 
-cdef struct func_state:
-    state st
-    random_uint64 f
+    cdef uint64_t splitmix64_next(splitmix64_state* state)  nogil
 
-ctypedef func_state func_state_t
 
+ctypedef uint64_t (*random_uint64)(splitmix64_state* state)
 
-cdef uint64_t _splitmix64(state_t *st):
-    cdef uint64_t z
-    # TODO: Use literals -- PyCharm complains
-    cdef uint64_t c1 = 11400714819323198485
-    cdef uint64_t c2 = 13787848793156543929
-    cdef uint64_t c3 = 10723151780598845931
-    st[0].state += c1 # 0x9E3779B97F4A7C15
-    z = <uint64_t>st[0].state
-    z = (z ^ (z >> 30)) * c2 # 0xBF58476D1CE4E5B9
-    z = (z ^ (z >> 27)) * c3 # 0x94D049BB133111EB
-    return z ^ (z >> 31)
+cdef uint64_t _splitmix64_anon(void* st) nogil:
+    return splitmix64_next(<splitmix64_state *>st)
 
-cdef uint64_t _splitmix64_anon(void* st):
-    return _splitmix64(<state *> st)
-
-cdef class CorePRNG:
+cdef class SplitMix64:
     """
     Prototype Core PRNG using directly implemented version of SplitMix64.
 
@@ -44,7 +30,7 @@ cdef class CorePRNG:
     Exposes no user-facing API except `get_state` and `set_state`. Designed
     for use in a `RandomGenerator` object.
     """
-    cdef state rng_state
+    cdef splitmix64_state rng_state
     cdef anon_func_state anon_func_state
     cdef public object _anon_func_state
 
@@ -70,7 +56,7 @@ cdef class CorePRNG:
         -----
         Testing only
         """
-        return _splitmix64(&self.rng_state)
+        return splitmix64_next(&self.rng_state)
 
     def get_state(self):
         """Get PRNG state"""

diff --git a/_randomgen/core_prng/src/splitmix64/splitmix64.c b/_randomgen/core_prng/src/splitmix64/splitmix64.c
@@ -0,0 +1,25 @@
+/*  Written in 2015 by Sebastiano Vigna (vigna@acm.org)
+
+To the extent possible under law, the author has dedicated all copyright
+and related and neighboring rights to this software to the public domain
+worldwide. This software is distributed without any warranty.
+
+See <http://creativecommons.org/publicdomain/zero/1.0/>.
+
+Modified 2018 by Kevin Sheppard.  Modifications licensed under the NCSA
+license.
+*/
+
+/* This is a fixed-increment version of Java 8's SplittableRandom generator
+   See http://dx.doi.org/10.1145/2714064.2660195 and
+   http://docs.oracle.com/javase/8/docs/api/java/util/SplittableRandom.html
+
+   It is a very fast generator passing BigCrush, and it can be useful if
+   for some reason you absolutely want 64 bits of state; otherwise, we
+   rather suggest to use a xoroshiro128+ (for moderately parallel
+   computations) or xorshift1024* (for massively parallel computations)
+   generator. */
+
+#include "splitmix64.h"
+
+extern inline uint64_t splitmix64_next(splitmix64_state *state);
diff --git a/_randomgen/core_prng/src/splitmix64/splitmix64.h b/_randomgen/core_prng/src/splitmix64/splitmix64.h
@@ -0,0 +1,12 @@
+#include <stdint.h>
+
+typedef struct s_splitmix64_state {
+  uint64_t state;
+} splitmix64_state;
+
+static inline uint64_t splitmix64_next(splitmix64_state *state) {
+  uint64_t z = (state->state += 0x9e3779b97f4a7c15);
+  z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
+  z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
+  return z ^ (z >> 31);
+}
diff --git a/_randomgen/core_prng/src/splitmix64/splitmix64.orig.c b/_randomgen/core_prng/src/splitmix64/splitmix64.orig.c
@@ -0,0 +1,28 @@
+/*  Written in 2015 by Sebastiano Vigna (vigna@acm.org)
+
+To the extent possible under law, the author has dedicated all copyright
+and related and neighboring rights to this software to the public domain
+worldwide. This software is distributed without any warranty.
+
+See <http://creativecommons.org/publicdomain/zero/1.0/>. */
+
+#include <stdint.h>
+
+/* This is a fixed-increment version of Java 8's SplittableRandom generator
+   See http://dx.doi.org/10.1145/2714064.2660195 and 
+   http://docs.oracle.com/javase/8/docs/api/java/util/SplittableRandom.html
+
+   It is a very fast generator passing BigCrush, and it can be useful if
+   for some reason you absolutely want 64 bits of state; otherwise, we
+   rather suggest to use a xoroshiro128+ (for moderately parallel
+   computations) or xorshift1024* (for massively parallel computations)
+   generator. */
+
+uint64_t x; /* The state can be seeded with any value. */
+
+uint64_t next() {
+	uint64_t z = (x += 0x9e3779b97f4a7c15);
+	z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
+	z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
+	return z ^ (z >> 31);
+}
diff --git a/_randomgen/core_prng/src/xoroshiro128/xoroshiro128.c b/_randomgen/core_prng/src/xoroshiro128/xoroshiro128.c
@@ -0,0 +1,53 @@
+/*  Written in 2016 by David Blackman and Sebastiano Vigna (vigna@acm.org)
+
+To the extent possible under law, the author has dedicated all copyright
+and related and neighboring rights to this software to the public domain
+worldwide. This software is distributed without any warranty.
+
+See <http://creativecommons.org/publicdomain/zero/1.0/>. */
+
+/* This is the successor to xorshift128+. It is the fastest full-period
+   generator passing BigCrush without systematic failures, but due to the
+   relatively short period it is acceptable only for applications with a
+   mild amount of parallelism; otherwise, use a xorshift1024* generator.
+
+   Beside passing BigCrush, this generator passes the PractRand test suite
+   up to (and included) 16TB, with the exception of binary rank tests, as
+   the lowest bit of this generator is an LFSR of degree 128. The next bit
+   can be described by an LFSR of degree 8256, but in the long run it will
+   fail linearity tests, too. The other bits needs a much higher degree to
+   be represented as LFSRs.
+
+   We suggest to use a sign test to extract a random Boolean value, and
+   right shifts to extract subsets of bits.
+
+   Note that the generator uses a simulated rotate operation, which most C
+   compilers will turn into a single instruction. In Java, you can use
+   Long.rotateLeft(). In languages that do not make low-level rotation
+   instructions accessible xorshift128+ could be faster.
+
+   The state must be seeded so that it is not everywhere zero. If you have
+   a 64-bit seed, we suggest to seed a splitmix64 generator and use its
+   output to fill s. */
+
+#include "xoroshiro128.h"
+
+extern inline uint64_t xoroshiro128_next(xoroshiro128_state *state);
+
+void xoroshiro128_jump(xoroshiro128_state *state) {
+  static const uint64_t JUMP[] = {0xbeac0467eba5facb, 0xd86b048b86aa9922};
+
+  uint64_t s0 = 0;
+  uint64_t s1 = 0;
+  for (int i = 0; i < sizeof JUMP / sizeof *JUMP; i++)
+    for (int b = 0; b < 64; b++) {
+      if (JUMP[i] & UINT64_C(1) << b) {
+        s0 ^= state->s[0];
+        s1 ^= state->s[1];
+      }
+      xoroshiro128_next(state);
+    }
+
+  state->s[0] = s0;
+  state->s[1] = s1;
+}
diff --git a/_randomgen/core_prng/src/xoroshiro128/xoroshiro128.h b/_randomgen/core_prng/src/xoroshiro128/xoroshiro128.h
@@ -0,0 +1,23 @@
+#include <stdint.h>
+
+typedef struct s_xoroshiro128_state {
+  uint64_t s[2];
+} xoroshiro128_state;
+
+static inline uint64_t rotl(const uint64_t x, int k) {
+  return (x << k) | (x >> (64 - k));
+}
+
+static inline uint64_t xoroshiro128_next(xoroshiro128_state *state) {
+  const uint64_t s0 = state->s[0];
+  uint64_t s1 = state->s[1];
+  const uint64_t result = s0 + s1;
+
+  s1 ^= s0;
+  state->s[0] = rotl(s0, 55) ^ s1 ^ (s1 << 14); // a, b
+  state->s[1] = rotl(s1, 36);                   // c
+
+  return result;
+}
+
+void xoroshiro128_jump(xoroshiro128_state *state);
diff --git a/_randomgen/core_prng/src/xoroshiro128/xoroshiro128plus.orig.c b/_randomgen/core_prng/src/xoroshiro128/xoroshiro128plus.orig.c
@@ -0,0 +1,74 @@
+/*  Written in 2016 by David Blackman and Sebastiano Vigna (vigna@acm.org)
+
+To the extent possible under law, the author has dedicated all copyright
+and related and neighboring rights to this software to the public domain
+worldwide. This software is distributed without any warranty.
+
+See <http://creativecommons.org/publicdomain/zero/1.0/>. */
+
+#include <stdint.h>
+
+/* This is the successor to xorshift128+. It is the fastest full-period
+   generator passing BigCrush without systematic failures, but due to the
+   relatively short period it is acceptable only for applications with a
+   mild amount of parallelism; otherwise, use a xorshift1024* generator.
+
+   Beside passing BigCrush, this generator passes the PractRand test suite
+   up to (and included) 16TB, with the exception of binary rank tests, as
+   the lowest bit of this generator is an LFSR of degree 128. The next bit
+   can be described by an LFSR of degree 8256, but in the long run it will
+   fail linearity tests, too. The other bits needs a much higher degree to
+   be represented as LFSRs.
+
+   We suggest to use a sign test to extract a random Boolean value, and
+   right shifts to extract subsets of bits.
+
+   Note that the generator uses a simulated rotate operation, which most C
+   compilers will turn into a single instruction. In Java, you can use
+   Long.rotateLeft(). In languages that do not make low-level rotation
+   instructions accessible xorshift128+ could be faster.
+
+   The state must be seeded so that it is not everywhere zero. If you have
+   a 64-bit seed, we suggest to seed a splitmix64 generator and use its
+   output to fill s. */
+
+uint64_t s[2];
+
+static inline uint64_t rotl(const uint64_t x, int k) {
+	return (x << k) | (x >> (64 - k));
+}
+
+uint64_t next(void) {
+	const uint64_t s0 = s[0];
+	uint64_t s1 = s[1];
+	const uint64_t result = s0 + s1;
+
+	s1 ^= s0;
+	s[0] = rotl(s0, 55) ^ s1 ^ (s1 << 14); // a, b
+	s[1] = rotl(s1, 36); // c
+
+	return result;
+}
+
+
+/* This is the jump function for the generator. It is equivalent
+   to 2^64 calls to next(); it can be used to generate 2^64
+   non-overlapping subsequences for parallel computations. */
+
+void jump(void) {
+	static const uint64_t JUMP[] = { 0xbeac0467eba5facb, 0xd86b048b86aa9922 };
+
+	uint64_t s0 = 0;
+	uint64_t s1 = 0;
+	for(int i = 0; i < sizeof JUMP / sizeof *JUMP; i++)
+		for(int b = 0; b < 64; b++) {
+			if (JUMP[i] & UINT64_C(1) << b) {
+				s0 ^= s[0];
+				s1 ^= s[1];
+			}
+			next();
+		}
+
+	s[0] = s0;
+	s[1] = s1;
+}