Merge pull request #20 from bashtage/dSFMT

ENH: Add initial support for dSFMT generator

Author: bashtage

README.md

@@ -9,6 +9,8 @@ The RNGs include:
 
 * [MT19937](https://github.com/numpy/numpy/blob/master/numpy/random/mtrand/),
  the NumPy rng
+* [dSFMT](http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/) a SSE2-aware 
+version of the MT19937 generator that is especially fast at generating doubles
 * [xorshift128+](http://xorshift.di.unimi.it/) and 
 [xorshift1024*](http://xorshift.di.unimi.it/)
 * [PCG32](http://www.pcg-random.org/) and [PCG64](http:w//www.pcg-random.org/)

randomstate/distributions.h

@@ -31,6 +31,8 @@ typedef int bool;
 #include "shims/mrg32k3a/mrg32k3a-shim.h"
 #elif defined(MLFG_1279_861_RNG)
 #include "shims/mlfg-1279-861/mlfg-1279-861-shim.h"
+#elif defined(DSFMT_RNG)
+#include "shims/dSFMT/dSFMT-shim.h"
 #else
 #error Unknown RNG!!!  Unknown RNG!!!  Unknown RNG!!!
 #endif

randomstate/interface.pyx

@@ -39,6 +39,8 @@ IF RNG_MOD_NAME == 'mrg32k3a':
     include "shims/mrg32k3a/mrg32k3a.pxi"
 IF RNG_MOD_NAME == 'mlfg_1279_861':
     include "shims/mlfg-1279-861/mlfg-1279-861.pxi"
+IF RNG_MOD_NAME == 'dsfmt':
+    include "shims/dSFMT/dSFMT.pxi"
 
 IF NORMAL_METHOD == 'inv':
     __normal_method = 'inv'

randomstate/performance.py

@@ -21,7 +21,7 @@
     scale_64 = 2
 
 RNGS = ['mlfg_1279_861', 'mrg32k3a', 'pcg64', 'pcg32', 'mt19937', 'xorshift128', 'xorshift1024',
-        'random']
+        'dsfmt', 'random']
 
 
 def timer(code, setup):

randomstate/prng/__init__.py

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

randomstate/prng/dsfmt/__init__.py

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

randomstate/shims/dSFMT/dSFMT-shim.c

@@ -0,0 +1,24 @@
+#include "dSFMT-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 void set_seed_by_array(aug_state* state, uint32_t init_key[], int key_length)
+{
+    dsfmt_init_by_array(state->rng, init_key, key_length);
+}
+
+void set_seed(aug_state* state, uint32_t seed)
+{
+    dsfmt_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]);
+}

randomstate/shims/dSFMT/dSFMT-shim.h

@@ -0,0 +1,54 @@
+#define RNG_TYPE rk_state
+
+#ifdef _WIN32
+#include "../../src/common/stdint.h"
+#define inline __inline
+#else
+#include <stdint.h>
+#endif
+
+#include "../../src/common/binomial.h"
+#include "../../src/common/entropy.h"
+#include "../../src/dSFMT/dSFMT.h"
+
+typedef struct s_aug_state {
+    dsfmt_t *rng;
+    binomial_t *binomial;
+
+    int has_gauss, shift_zig_random_int, has_uint32;
+    double gauss;
+    uint32_t uinteger;
+    uint64_t zig_random_int;
+} aug_state;
+
+static inline uint32_t random_uint32(aug_state* state)
+{
+    double d = dsfmt_genrand_close1_open2(state->rng);
+    uint64_t *out = (uint64_t *)&d;
+    return (uint32_t)(*out & 0xffffffff);
+}
+
+static inline uint64_t random_uint64(aug_state* state)
+{
+    double d = dsfmt_genrand_close1_open2(state->rng);
+    uint64_t out;
+    uint64_t *tmp;
+    tmp = (uint64_t *)&d;
+    out = *tmp << 32;
+    d = dsfmt_genrand_close1_open2(state->rng);
+    tmp =  (uint64_t *)&d;
+    out |= *tmp & 0xffffffff;
+    return out;
+}
+
+static inline double random_double(aug_state* state)
+{
+    return dsfmt_genrand_close1_open2(state->rng) - 1.0;
+}
+
+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);
+

randomstate/shims/dSFMT/dSFMT.pxi

@@ -0,0 +1,104 @@
+DEF RNG_NAME = 'dSFMT'
+DEF RNG_ADVANCEABLE = 0
+DEF RNG_SEED = 1
+DEF RNG_JUMPABLE = 0
+DEF RNG_STATE_LEN = 4
+DEF NORMAL_METHOD = 'zig'
+DEF DSFMT_MEXP = 19937
+DEF DSFMT_N = 191 # ((DSFMT_MEXP - 128) / 104 + 1)
+DEF DSFMT_N_PLUS_1 = 192 # DSFMT_N + 1
+
+ctypedef uint32_t rng_state_t
+
+cdef extern from "distributions.h":
+
+    cdef union W128_T:
+        uint64_t u[2];
+        uint32_t u32[4];
+        double d[2];
+
+    ctypedef W128_T w128_t;
+
+    cdef struct DSFMT_T:
+        w128_t status[DSFMT_N_PLUS_1];
+        int idx;
+
+    ctypedef DSFMT_T dsfmt_t;
+
+    cdef struct s_aug_state:
+        dsfmt_t *rng
+        binomial_t *binomial
+
+        int has_gauss, shift_zig_random_int, has_uint32
+        double gauss
+        uint64_t zig_random_int
+        uint32_t uinteger
+
+    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)
+
+ctypedef dsfmt_t rng_t
+
+cdef object _get_state(aug_state state):
+    cdef uint32_t [:] key = np.zeros(4 * DSFMT_N_PLUS_1, dtype=np.uint32)
+    cdef Py_ssize_t i, j, key_loc = 0
+    cdef w128_t state_val
+    for i in range(DSFMT_N_PLUS_1):
+        state_val = state.rng.status[i]
+        for j in range(4):
+            key[key_loc] = state_val.u32[j]
+            key_loc += 1
+    return (np.asarray(key), state.rng.idx)
+
+cdef object _set_state(aug_state state, object state_info):
+    cdef uint32_t [:] key = state_info[0]
+    cdef Py_ssize_t i, j, key_loc = 0
+    for i in range(DSFMT_N_PLUS_1):
+        for j in range(4):
+            state.rng.status[i].u32[j] = key[key_loc]
+            key_loc += 1
+
+    state.rng.idx = state_info[1]
+
+DEF CLASS_DOCSTRING = """
+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. 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}, optional
+    Random seed initializing the pseudo-random number generator.
+    Can be an integer in [0, 2**32] or ``None`` (the default).
+    If `seed` is ``None``, then `dSFMT.RandomState` will try to read
+    entropy from ``/dev/urandom`` (or the Windows analogue) 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`. `RandomState`, besides being
+NumPy-aware, has the advantage that it provides a much larger number
+of probability distributions to choose from.
+"""

randomstate/src/dSFMT/LICENSE.txt

@@ -0,0 +1,32 @@
+Copyright (c) 2007, 2008, 2009 Mutsuo Saito, Makoto Matsumoto
+and Hiroshima University.
+Copyright (c) 2011, 2002 Mutsuo Saito, Makoto Matsumoto, Hiroshima
+University and The University of Tokyo.
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+    * Redistributions in binary form must reproduce the above
+      copyright notice, this list of conditions and the following
+      disclaimer in the documentation and/or other materials provided
+      with the distribution.
+    * Neither the name of the Hiroshima University nor the names of
+      its contributors may be used to endorse or promote products
+      derived from this software without specific prior written
+      permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.