Merge pull request #97 from bashtage/add-sfmt

ENH: Add SFMT PRNG

Author: bashtage

randomstate/distributions.h

@@ -40,6 +40,8 @@ typedef int bool;
 #include "interface/mlfg-1279-861/mlfg-1279-861-shim.h"
 #elif defined(RS_DSFMT)
 #include "interface/dSFMT/dSFMT-shim.h"
+#elif defined(RS_SFMT)
+#include "interface/sfmt/sfmt-shim.h"
 #else
 #error Unknown RNG!!!  Unknown RNG!!!  Unknown RNG!!!
 #endif

randomstate/interface/pcg-64/pcg-64-shim.h

@@ -65,7 +65,5 @@ inline void entropy_init(aug_state* state)
 
 inline double random_double(aug_state* state)
 {
-    uint64_t rn;
-    rn = random_uint64(state);
-    return (rn >> 11) * (1.0 / 9007199254740992.0);
+    return (random_uint64(state) >> 11) * (1.0 / 9007199254740992.0);
 }

randomstate/interface/sfmt/sfmt-shim.c

@@ -0,0 +1,39 @@
+#include "sfmt-shim.h"
+
+extern inline uint32_t random_uint32(aug_state *state);
+
+extern inline uint64_t random_uint64(aug_state *state);
+
+extern inline double random_double(aug_state *state);
+
+extern inline uint64_t random_raw_values(aug_state *state);
+
+void reset_buffer(aug_state *state) {
+  int i = 0;
+  for (i = 0; i < (2 * SFMT_N); i++) {
+    state->buffered_uint64[i] = 0ULL;
+  }
+  state->buffer_loc = 2 * SFMT_N;
+}
+
+extern void set_seed_by_array(aug_state *state, uint32_t init_key[],
+                              int key_length) {
+  reset_buffer(state);
+  sfmt_init_by_array(state->rng, init_key, key_length);
+}
+
+void set_seed(aug_state *state, uint32_t seed) {
+  reset_buffer(state);
+  sfmt_init_gen_rand(state->rng, seed);
+}
+
+void entropy_init(aug_state *state) {
+  uint32_t seeds[1];
+  entropy_fill((void *)seeds, sizeof(seeds));
+  set_seed(state, seeds[0]);
+}
+
+// void jump_state(aug_state* state)
+//{
+//    dSFMT_jump(state->rng, poly_128);
+//}

randomstate/interface/sfmt/sfmt-shim.h

@@ -0,0 +1,68 @@
+#ifdef _WIN32
+#include "../../src/common/stdint.h"
+#define inline __forceinline
+#else
+#include <stdint.h>
+#endif
+
+#include "../../src/common/binomial.h"
+#include "../../src/common/entropy.h"
+#include "../../src/sfmt/sfmt.h"
+
+typedef struct s_aug_state {
+  sfmt_t *rng;
+  binomial_t *binomial;
+
+  int has_gauss, has_gauss_float, shift_zig_random_int, has_uint32;
+  float gauss_float;
+  double gauss;
+  uint32_t uinteger;
+  uint64_t zig_random_int;
+
+  uint64_t *buffered_uint64;
+  int buffer_loc;
+} aug_state;
+
+static inline uint64_t random_uint64_from_buffer(aug_state *state) {
+  uint64_t out;
+  if (state->buffer_loc >= (2 * SFMT_N)) {
+    state->buffer_loc = 0;
+    sfmt_fill_array64(state->rng, state->buffered_uint64, 2 * SFMT_N);
+  }
+  out = state->buffered_uint64[state->buffer_loc];
+  state->buffer_loc++;
+  return out;
+}
+
+static inline uint32_t random_uint32(aug_state *state) {
+  uint64_t d;
+  if (state->has_uint32) {
+    state->has_uint32 = 0;
+    return state->uinteger;
+  }
+  d = random_uint64_from_buffer(state);
+  state->uinteger = (uint32_t)(d >> 32);
+  state->has_uint32 = 1;
+  return (uint32_t)(d & 0xFFFFFFFFUL);
+}
+
+static inline uint64_t random_uint64(aug_state *state) {
+  return random_uint64_from_buffer(state);
+}
+
+static inline double random_double(aug_state *state) {
+  return (random_uint64_from_buffer(state) >> 11) * (1.0 / 9007199254740992.0);
+}
+
+static inline uint64_t random_raw_values(aug_state *state) {
+  return random_uint64_from_buffer(state);
+}
+
+extern void entropy_init(aug_state *state);
+
+extern void set_seed_by_array(aug_state *state, uint32_t init_key[],
+                              int key_length);
+
+extern void set_seed(aug_state *state, uint32_t seed);
+
+// extern void jump_state(aug_state* state);

randomstate/interface/sfmt/sfmt.pxi

@@ -0,0 +1,142 @@
+DEF RS_RNG_NAME = u'sfmt'
+DEF SFMT_MEXP = 19937
+DEF SFMT_N = 156 # (SFMT_MEXP / 128 + 1)
+DEF RS_SEED_NBYTES = 1
+DEF RS_SEED_ARRAY_BITS = 32
+
+ctypedef uint32_t rng_state_t
+
+cdef extern from "distributions.h":
+
+    cdef union W128_T:
+        uint32_t u[4]
+        uint64_t u64[2]
+
+    ctypedef W128_T w128_t;
+
+    cdef struct SFMT_T:
+        w128_t state[SFMT_N];
+        int idx;
+
+    ctypedef SFMT_T sfmt_t;
+
+    cdef struct s_aug_state:
+        sfmt_t *rng
+        binomial_t *binomial
+
+        int has_gauss, shift_zig_random_int, has_uint32, has_gauss_float
+        float gauss_float
+        double gauss
+        uint64_t zig_random_int
+        uint32_t uinteger
+
+        uint64_t *buffered_uint64
+        int buffer_loc
+
+    ctypedef s_aug_state aug_state
+
+    cdef void set_seed(aug_state* state, uint32_t seed)
+
+    cdef void set_seed_by_array(aug_state* state, uint32_t init_key[], int key_length)
+
+    cdef void jump_state(aug_state* state)
+
+
+ctypedef sfmt_t rng_t
+
+cdef object _get_state(aug_state state):
+    cdef uint32_t [::1] key = np.zeros(4 * SFMT_N, dtype=np.uint32)
+    cdef uint64_t [::1] buf = np.zeros(2 * SFMT_N, dtype=np.uint64)
+    cdef Py_ssize_t i, j, key_loc = 0
+    cdef w128_t state_val
+    for i in range(SFMT_N):
+        state_val = state.rng.state[i]
+        for j in range(4):
+            key[key_loc] = state_val.u[j]
+            key_loc += 1
+    for i in range(2 * SFMT_N):
+        buf[i] = state.buffered_uint64[i]
+
+    return (np.asarray(key), state.rng.idx,
+            np.asarray(buf), state.buffer_loc)
+
+cdef object _set_state(aug_state *state, object state_info):
+    cdef Py_ssize_t i, j, key_loc = 0
+    cdef uint32_t [::1] key = state_info[0]
+    state.rng.idx = state_info[1]
+
+    for i in range(SFMT_N):
+        for j in range(4):
+            state.rng.state[i].u[j] = key[key_loc]
+            key_loc += 1
+
+    state.buffer_loc = <int>state_info[3]
+    for i in range(2 * SFMT_N):
+        state.buffered_uint64[i] = state_info[2][i]
+
+
+DEF CLASS_DOCSTRING = u"""
+RandomState(seed=None)
+
+Container for the SIMD-based Mersenne Twister pseudo-random number generator.
+
+``dSFMT.RandomState`` exposes a number of methods for generating random
+numbers drawn from a variety of probability distributions [1]_ . In addition
+to the distribution-specific arguments, each method takes a keyword argument
+`size` that defaults to ``None``. If `size` is ``None``, then a single
+value is generated and returned. If `size` is an integer, then a 1-D
+array filled with generated values is returned. If `size` is a tuple,
+then an array with that shape is filled and returned.
+
+**No Compatibility Guarantee**
+
+``dSFMT.RandomState`` does not make a guarantee that a fixed seed and a
+fixed series of calls to ``dSFMT.RandomState`` methods using the same
+parameters will always produce the same results. This is different from
+``numpy.random.RandomState`` guarantee. This is done to simplify improving
+random number generators.  To ensure identical results, you must use the
+same release version.
+
+Parameters
+----------
+seed : {None, int, array_like}, optional
+    Random seed initializing the pseudo-random number generator.
+    Can be an integer in [0, 2**32-1], array of integers in
+    [0, 2**32-1] or ``None`` (the default). If `seed` is ``None``,
+    then ``dSFMT.RandomState`` will try to read entropy from
+    ``/dev/urandom`` (or the Windows analog) if available to
+    produce a 64-bit seed. If unavailable, the a 64-bit hash of the time
+    and process ID is used.
+
+Notes
+-----
+The Python stdlib module "random" also contains a Mersenne Twister
+pseudo-random number generator with a number of methods that are similar
+to the ones available in ``RandomState``. The `RandomState` object, besides
+being NumPy-aware, also has the advantage that it provides a much larger
+number of probability distributions to choose from.
+
+>>> from randomstate.entropy import random_entropy
+>>> import randomstate.prng.sfmt as rnd
+>>> seed = random_entropy()
+>>> rs = [rnd.RandomState(seed) for _ in range(10)]
+
+**State and Seeding**
+
+The ``sfmt.RandomState`` state vector consists of a 624 element array of
+32-bit unsigned integers plus a single integer value between 0 and 624
+indicating  the current position within the main array. The implementation
+used here augments this with a array of 64-bit unsigned integers which are used
+to efficiently access the random numbers produced by the SFMT generator.
+
+``sfmt.RandomState`` is seeded using either a single 32-bit unsigned integer
+or a vector of 32-bit unsigned integers.  In either case, the input seed is
+used as an input (or inputs) for a hashing function, and the output of the
+hashing function is used as the initial state. Using a single 32-bit value
+for the seed can only initialize a small range of the possible initial
+state values.
+
+.. [1] Mutsuo Saito and Makoto Matsumoto, "SIMD-oriented Fast Mersenne
+       Twister: a 128-bit Pseudorandom Number Generator." Monte Carlo
+       and Quasi-Monte Carlo Methods 2006, Springer, pp. 607 -- 622, 2008.
+"""

randomstate/interface/xorshift1024/xorshift1024-shim.h

@@ -47,9 +47,7 @@ inline uint64_t random_raw_values(aug_state* state)
 
 inline double random_double(aug_state* state)
 {
-    uint64_t rn;
-    rn = random_uint64(state);
-    return (rn >> 11) * (1.0 / 9007199254740992.0);
+    return (random_uint64(state) >> 11) * (1.0 / 9007199254740992.0);
 }
 
 extern void set_seed(aug_state* state, uint64_t seed);

randomstate/interface/xorshift128/xorshift128-shim.h

@@ -47,9 +47,7 @@ inline uint64_t random_raw_values(aug_state* state)
 
 inline double random_double(aug_state* state)
 {
-    uint64_t rn;
-    rn = random_uint64(state);
-    return (rn >> 11) * (1.0 / 9007199254740992.0);
+    return (random_uint64(state) >> 11) * (1.0 / 9007199254740992.0);
 }
 
 extern void set_seed(aug_state* state, uint64_t seed);

randomstate/performance.py

@@ -23,7 +23,7 @@
     scale_64 = 2
 
 RNGS = ['mlfg_1279_861', 'mrg32k3a', 'pcg64', 'pcg32', 'mt19937', 'xorshift128', 'xorshift1024',
-        'xoroshiro128plus', 'dsfmt', 'random']
+        'xoroshiro128plus', 'dsfmt', 'sfmt', 'random']
 
 
 def timer(code, setup):
@@ -41,6 +41,7 @@ def run_timer(dist, command, numpy_command=None, setup='', random_type=''):
 
     res = {}
     for rng in RNGS:
+        print(str(rng))
         mod = 'randomstate.prng' if rng != 'random' else 'numpy'
         key = '-'.join((mod, rng, dist)).replace('"', '')
         command = numpy_command if 'numpy' in mod else command

randomstate/prng/__init__.py

@@ -7,6 +7,7 @@
 from .pcg32 import pcg32
 from .pcg64 import pcg64
 from .dsfmt import dsfmt
+from .sfmt import sfmt
 
 def __generic_ctor(mod_name='mt19937'):
     """
@@ -46,6 +47,8 @@ def __generic_ctor(mod_name='mt19937'):
         mod = xorshift1024
     elif mod_name == 'dsfmt':
         mod = dsfmt
+    elif mod_name == 'sfmt':
+        mod = sfmt
     else:
         raise ValueError(str(mod_name) + ' is not a known PRNG module.')
 

randomstate/prng/sfmt/__init__.py

@@ -0,0 +1 @@
+from .sfmt import *