ENH: Add jumt to MT19937

Add jump function to MT19937 generator using Horner method

closes #84

Author: bashtage

.travis.yml

@@ -21,7 +21,7 @@ matrix:
   - env:
     - PYTHON=3.4
     - NUMPY=1.10
-    - CYTHON=0.23
+    - CYTHON=0.24
   - env:
     - PYTHON=3.5
   - env:

README.md

@@ -144,7 +144,7 @@ Building requires:
 
   * Python (2.7, 3.4, 3.5, 3.6)
   * NumPy (1.9, 1.10, 1.11, 1.12)
-  * Cython (0.22, 0.23, 0.24, 0.25)
+  * Cython (0.22, **not** 0.23, 0.24, 0.25)
   * tempita (0.5+), if not provided by Cython
 
 Testing requires pytest (3.0+).

README.rst

@@ -160,7 +160,7 @@ Building requires:
 
 -  Python (2.7, 3.4, 3.5, 3.6)
 -  NumPy (1.9, 1.10, 1.11, 1.12)
--  Cython (0.22, 0.23, 0.24, 0.25)
+-  Cython (0.22, **not** 0.23, 0.24, 0.25)
 -  tempita (0.5+), if not provided by Cython
 
 Testing requires pytest (3.0+).

doc/source/change-log.rst

@@ -33,6 +33,7 @@ Since Version 1.13
 * Add SIMD-oriented Fast Mersenne Twister
   (`SFMT <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/>`_) generator.
 * Add complex normal (:meth:`~randomstate.prng.mt19937.complex_normal`)
+* Added support for jumping the MT19937 generator
 
 Version 1.13
 ------------

doc/source/index.rst

@@ -88,7 +88,9 @@ The main innovation is the inclusion of a number of alternative pseudo-random nu
 generators, 'in addition' to the standard PRNG in NumPy.  The included PRNGs are:
 
 * MT19937 - The standard NumPy generator.  Produces identical results to NumPy
-  using the same seed/state. See `NumPy's documentation`_.
+  using the same seed/state. Adds a jump function that advances the generator
+  as-if 2**128 draws have been made (:meth:`randomstate.prng.mt19937.jump`).
+  See `NumPy's documentation`_.
 * dSFMT - A SSE2 enables version of the MT19937 generator.  Theoretically the
   same, but with a different state and so it is not possible to produce a
   sequence identical to MT19937. See the `dSFMT authors' page`_.

doc/source/mt19937.rst

@@ -16,6 +16,12 @@ Random generator
    get_state
    set_state
 
+Parallel generation
+===================
+.. autosummary::
+   :toctree: generated/
+
+   jump
 
 Simple random data
 ==================

doc/source/parallel.rst

@@ -54,6 +54,8 @@ are listed below.
 +-----------------+-------------------------+-------------------------+-------------------------+
 | xorshift1024    | :math:`2^{1024}`        | :math:`2^{512}`         | 64                      |
 +-----------------+-------------------------+-------------------------+-------------------------+
+| mt19937         | :math:`2^{19937}`        | :math:`2^{128}`        | 32                      |
++-----------------+-------------------------+-------------------------+-------------------------+
 
 ``jump`` can be used to produce long blocks which should be long enough to not
 overlap.

pyproject.toml

@@ -0,0 +1,2 @@
+[build-system]
+requires = ["numpy>=1.09", "cython>=0.24", "setuptools", "wheel"]

randomstate/interface/random-kit/random-kit-poly.h

@@ -1,4 +1,6 @@
 #include "random-kit-shim.h"
+#include "random-kit-poly.h"
+
 
 extern inline uint32_t random_uint32(aug_state* state);
 
@@ -23,4 +25,9 @@ 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)
+{
+    randomkit_jump(state->rng, poly);
 }