ENH: jump added to dSFMT

Sheppard, Kevin · bashtage · commit 1078941ea884 · 2016-01-19T21:15:35.000-05:00
Add jump to dSFMT that jumps the rng 2**128 moves
diff --git a/doc/source/dsfmt.rst b/doc/source/dsfmt.rst
@@ -15,6 +15,12 @@ Random generator
    get_state
    set_state
 
+Parallel generation
+===================
+.. autosummary::
+   :toctree: generated/
+
+   jump
 
 Simple random data
 ==================
diff --git a/randomstate/shims/dSFMT/dSFMT-poly.h b/randomstate/shims/dSFMT/dSFMT-poly.h
@@ -0,0 +1,53 @@
+static const char * poly_128 =
+"f4dfa6c62049d0776e0bf6f1e953f3aa38abb113df86be024eab3773ad5f2b82ead936022e656dff7e562691c59dd5f7d2"
+"566b78d9669002503c4ddb1888a49f32333f515e6c60c4ecd221078ec6f26f0a90f4875067ca1f399a99775037adf90556"
+"6e2c7e6b42131420f8f04f112c92621c9b1502f2a8aefad6c667904af62f0d55e02d396902d3b89450103c5ce5fe0408d9"
+"7cbb864861b49e4e42048ff3310b48faac55095a7f422eea4aade752f947f947c6be0a0c665bdea099246ab9eff658ea8c"
+"a468bf49d0227748367878de06d7bd86ea6708fcac6e252f5f00f04309b2aac3036b64afb39d990427c6c9f03477cc7e93"
+"5c43c0e61bc161db8eb15516eee8cb377ecbc1849207990fb6778721b29bfe0d89bfda1b3772fa5b0b1f7ec3daf3605203"
+"2285898c6f6396f55010c31f8201b7e2e51d94f920bfe57684c5415cc342cb39a0045d9793d13cf8646096daeb8bb9bfc2"
+"0a90de8f2426da8733267a9b9674f32154e8f84a9932223a2ca3c787d0b66df6675febbdfcba2f9cef09c621c57e11098b"
+"3289c77397aaae8b104642ffe0c4b75598efbc53745984d68b4d6656cae299ae2be55217a9a02b009ca7be32f47fbe434b"
+"ce4914a34d0c9b0085bede9b8a99319c34660d66f0124b5a7714c4bf3cbfec3ee43ed817087168bad80133bebaeeb68cf7"
+"929a24d1bb3de831a8340d220906ab04159cf94b21d5ee813bd7c80f10f01b43052af530917513b169254c25d6fcfe6cb4"
+"20d6ce92f54886ef6eaf9a5ba35e893ff593834d05ddf28899e42d729c7df3d21ef036020789739366f0c11ec52ff92a0b"
+"fd8ba69508e27b20fabb8217bd36b90e5aa918159ac87913bc7b46c04e366c23c92807fbe9c6a407e6a4db0b4fc23c3b6c"
+"706b5ca058fe8c190f849f18d16d6b48b5ed760eb202fd566291a799420b9654e08b8118bcbfead8e9dd2fdb9b053e9bdf"
+"b665285c78718f726d0b3d6c37e116428ec9ac9db2637259e4e8d6402bbada46c6bdb03985e19a82e9b4e57de1b025a3cb"
+"1f850beae7e8da9941655825bce0e89d536b6ee9064865b1a85c185e9fc9cb7f435de13d44773c00eed442a286e4ab807e"
+"3cab4dc3441d1b7d2af693812ae8b39652bb8c835fc895d13d6da93541afeadeee450475c29f3b2dfa8ef1c1e2547463b2"
+"cc2f0ff7a42ac4dd35e25c4fa030d2d2766fbe9f2d04c1304671747bace2f7dd55142bfa60f8cbc968bfc3d7a342152dc6"
+"84a0fb5a32c0962a62b5220ac0f72add9d8b84d6cc76b97d03245e01fc8da3414a49bb4075d3488f29b56dc42ba69e3b58"
+"529448c943ecfd98b3784a39d0b8609a8fb945e757f4569f53bd2cf80f7f638acf5b67fe9c560a3b7b0cf7e0398f31aa8b"
+"03cf9c62b24296b6d8596b694469a02686c38daa16a1ef86e012d61a2f7de1693a5c00b3685175caec3c67146477eba548"
+"30f1d546cb18a553779aa46adb4f2010e33f3def847c7d89b51a8462b227605f6c920fd558a6daf64bc98682e508ae960c"
+"0c571870e603ba1fce0c13d53176f353fd319959e13db93eae1359f06e3dd4767c04f824cf34ec7bf8f60161ba1a615db8"
+"2852eca9e3869afa711ab9a090660b0dc6cfbea310dda77e02310fbaeacd2636f975838c2dbcdbe9ac2cd85cee28f5e3f0"
+"c73abf62f9fa02cd79a7606b7ba855db68a07848b057c3aaf38f1a70086e14616f6f88305a1f9ce6b41378a620d4db3e0e"
+"7e1d421590dccaeff86212e232eeb5eb8a8d33a8c9b25ae88f3a7bd5032b4efa68f8af3186a02ffcbf5456f12beccace94"
+"c81c360cc4a0dcc642b59f991eec68c59af78139ca60b96d6a18e9535f8995e89bd2cf6a0aef3acffd33d1c0c1b79b6641"
+"4a91d9f65b2b4ec65844b96f725d2b4b0c309f3eb9d714e9dd939bbdfd85ce8fb43679aeab13f6c29549949503c9466dbd"
+"337c4cdde46d6eacd15f21f4d8fdeaa627a47884c88a9c85f0b731d271a8ea7cb9e04a4a149c23c10f56b3a0476dc77a99"
+"9d6e4f813e4b0f805e2a693e2ae4ae0ecc423c9ba5d17b42e691abf83784a582f2b1fd85d1e0a27ba38a500963568b2450"
+"363d2c5e3f7b8ba3e5b56e4e9f745a3a710bf2ae233c303068c532ce78ff031e6ab28b705dd94d7db4500909edb5626b8c"
+"9bd5ff4f0b4741388f0b91563ee516934c013e901572cba005ac5c535f4f107903be9af7b2793dfb61b5070facbe71eefe"
+"1b5600f975c8c38c3a2350d78beadfecb78e981164ae8bc866e732972d3ceef4aac68e15861f9b881d9b51b4edece150bc"
+"124b07645defb4202ef5d0e0962db98cae6ed459561c93c74c20bd64362e4f4fffc389a6cd80514604ff22eecc10c9cbc7"
+"981d19a8102b24146354c463107c9dc070e29e70df3578022acf72289ef071ab9f9402a544d0399f1b1e5f206b6d46d445"
+"f6d612a490e72918e00c853eda8493bef511149e80c9ab56e8b4b8cba3987249f77d060e61760e5792ac321c987c03c260"
+"6e9393a7970212992cdbd16448078d5039d4c2c3199714f53278f4f7b1d2e514cf95bdfc078b8bb0db659cb2c3f5cc0289"
+"0ea84f05d414c88d2db9e9f8455659b9fa6254405317245fa070d6970cafb4dadb2522b490a5c8e02fe973a8cdbfbfbdbf"
+"b01535099ffba3d3896bc4d1189fc570c3e6fdc6469265b8da912772e75dd62ab71be507f700d56cac5e68fd6b57ec1661"
+"68ab5258a69625c142a5b1b3519f94be1bde5e51d3bd8ea0c12d5af2fe4615b1b7bd4a96628a4fabc65925ff09718f63bb"
+"ebaad98f89bd9543a27b3ff3b5d8bfa89f941a5eb8cc005ccd4a705190e1c9dc6a9f4264e5ee658520a4438e92de854bff"
+"c39f8dc7dfbb5de4f14ba63ea16a37d14a7b4610f95b6cffd55e4679b29cedbdf20e7bd16da822fad910c359ee3a68e48a"
+"ae6e769b0e291d5d3aa3e2ca9d8d23abe8a1d5349f4991e9300852cc0befb20c2fc0d169306b260763344024f8092cbcc2"
+"4c6807363e9fc548a30d5faab3a94b2af0782a2942be80c45d8b0587efd587394ef33c33022436e285806ddffdd32fe363"
+"45c3c38ed8d680abeb7a028b44ee6f94d060a14c7019bb6af1f1b5f0a562957d19826d8cc216f9b908c989ccd5415e3525"
+"dfe9422ffb5b50b7cc3083dc325544751e5683535d7439d3da2b0bb73bea551dd99e04e0e793804f4774eb6b1daf781d9c"
+"aa5128274e599e847862fe309027813d3e4eda0bbeb7201856a5c5d8370e44dabff0bb229c723ba0a6bcf29c44536147de"
+"11b7835991018100105bd4329217f7386903fe8e7363cd7b3e893244e245e0a187467664c05b0be1fd429722b9b9a5e319"
+"8147fad72776e8a63aab9054fa9d259af0198d088d71d132e6068676a8e9ebb0f616b51ee34aac39c2c2221c7112401727"
+"0d75ff4a048363c389e04e9b440ad2032a381ac2cfc54f409caa791e65ee4f5d6cd035008f219b88a803a7382ae447bf65"
+"a3df2176b25b3b7b67dabe34decd9a1384dc7a003916ca8fbcb29b3ad6fd8eac5bbbaa3bdfa6c6a3ad9427c4f3ed79fea2"
+"6e14c8ce5fa3b4f82c5f7b6d2125916753a7b92ce9b46d45";
diff --git a/randomstate/shims/dSFMT/dSFMT-shim.c b/randomstate/shims/dSFMT/dSFMT-shim.c
@@ -1,4 +1,5 @@
 #include "dSFMT-shim.h"
+#include "dSFMT-poly.h"
 
 extern inline uint32_t random_uint32(aug_state* state);
 
@@ -33,4 +34,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)
+{
+    dSFMT_jump(state->rng, poly_128);
 }
diff --git a/randomstate/shims/dSFMT/dSFMT-shim.h b/randomstate/shims/dSFMT/dSFMT-shim.h
@@ -9,8 +9,10 @@
 
 #include "../../src/common/binomial.h"
 #include "../../src/common/entropy.h"
+#include "../../src/dSFMT/dSFMT-jump.h"
 #include "../../src/dSFMT/dSFMT.h"
 
+
 typedef struct s_aug_state {
     dsfmt_t *rng;
     binomial_t *binomial;
@@ -69,3 +71,4 @@ extern void set_seed_by_array(aug_state* state, uint32_t init_key[], int key_len
 
 extern void set_seed(aug_state* state, uint32_t seed);
 
+extern void jump_state(aug_state* state);
diff --git a/randomstate/shims/dSFMT/dSFMT.pxi b/randomstate/shims/dSFMT/dSFMT.pxi
@@ -1,4 +1,5 @@
 DEF RS_RNG_NAME = 'dSFMT'
+DEF RS_RNG_JUMPABLE = 1
 DEF DSFMT_MEXP = 19937
 DEF DSFMT_N = 191 # ((DSFMT_MEXP - 128) / 104 + 1)
 DEF DSFMT_N_PLUS_1 = 192 # DSFMT_N + 1
@@ -38,6 +39,9 @@ cdef extern from "distributions.h":
 
     cdef void set_seed_by_array(aug_state* state, uint32_t init_key[], int key_length)
 
+    cdef void jump_state(aug_state* state)
+
+
 ctypedef dsfmt_t rng_t
 
 cdef object _get_state(aug_state state):
@@ -80,8 +84,8 @@ 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
+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,
@@ -113,4 +117,28 @@ 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.
+
+**Parallel Features**
+
+``dsfmt.RandomState`` can be used in parallel applications by
+calling the method ``jump`` which advances the the state as-if :math:`2^{128}`
+random numbers have been generated [2]_. This allow the original sequence to
+be split so that distinct segments can be used on each worker process.  All
+generators should be initialized with the same seed to ensure that the
+segments come from the same sequence.
+
+>>> from randomstate.entropy import random_entropy
+>>> import randomstate.prng.dsfmt as rnd
+>>> seed = random_entropy()
+>>> rs = [rnd.RandomState(seed) for _ in range(10)]
+# Advance rs[i] by i jumps
+>>> for i in range(10):
+        rs[i].jump(i)
+
+.. [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.
+.. [2] Hiroshi Haramoto, Makoto Matsumoto, and Pierre L'Ecuyer, "A Fast
+       Jump Ahead Algorithm for Linear Recurrences in a Polynomial Space",
+       Sequences and Their Applications - SETA, 290--298, 2008.
 """
diff --git a/randomstate/shims/xorshift1024/xorshift1024.pxi b/randomstate/shims/xorshift1024/xorshift1024.pxi
@@ -50,7 +50,7 @@ RandomState(seed=None)
 Container for the xorshift1024* pseudo random number generator.
 
 xorshift1024* is a 64-bit implementation of Saito and Matsumoto's XSadd
-generator [1]_,[2]_,[3]_,[4]_. xorshift1024* has a period of 2**1024 - 1 and
+generator [1]_. xorshift1024* has a period of 2**1024 - 1 and
 supports jumping the sequence in increments of 2**512, which allow multiple
 non-overlapping sequences to be generated.
 
@@ -90,7 +90,8 @@ See xorshift128 for a faster implementation that has a smaller period.
 calling the method ``jump`` which advances the
 the state as-if :math:`2^{512}` random numbers have been generated. This
 allow the original sequence to be split so that distinct segments can be used
-on each worker process.
+on each worker process. All generators should be initialized with the same
+seed to ensure that the segments come from the same sequence.
 
 >>> import randomstate.prng.xorshift1024 as rnd
 >>> rs = [rnd.RandomState(1234) for _ in range(10)]
diff --git a/randomstate/shims/xorshift128/xorshift128.pxi b/randomstate/shims/xorshift128/xorshift128.pxi
@@ -45,7 +45,7 @@ RandomState(seed=None)
 Container for the xorshift128+ pseudo random number generator.
 
 xorshift128+ is a 64-bit implementation of Saito and Matsumoto's XSadd
-generator [1]_,[2]_,[3]_,[4]_. xorshift128+ has a period of 2**128 - 1 and
+generator [1]_. xorshift128+ has a period of 2**128 - 1 and
 supports jumping the sequence in increments of 2**64, which allow multiple
 non-overlapping sequences to be generated.
 
@@ -86,7 +86,8 @@ See xorshift1024 for an implementation with a larger period
 calling the method ``jump`` which advances the
 the state as-if :math:`2^{64}` random numbers have been generated. This
 allow the original sequence to be split so that distinct segments can be used
-on each worker process.
+on each worker process. All generators should be initialized with the same
+seed to ensure that the segments come from the same sequence.
 
 >>> import randomstate.prng.xorshift128 as rnd
 >>> rs = [rnd.RandomState(1234) for _ in range(10)]
diff --git a/randomstate/src/dSFMT/dSFMT-jump.c b/randomstate/src/dSFMT/dSFMT-jump.c
@@ -49,7 +49,7 @@ extern "C" {
 		  UINT64_C(0x3ffafa9bfbd5d55d),
 		  UINT64_C(0x0470470470573573),
 		  UINT64_C(0x0250250251259259)}},
-	{0, {}}
+	{0}
     };
 
     inline static void next_state(dsfmt_t * dsfmt);
@@ -150,11 +150,13 @@ extern "C" {
 	dsfmt_t work;
 	int index = dsfmt->idx;
 	int bits;
+	int i;
+	int j;
 	memset(&work, 0, sizeof(dsfmt_t));
 	add_fix(dsfmt);
 	dsfmt->idx = DSFMT_N64;
 
-	for (int i = 0; jump_string[i] != '\0'; i++) {
+	for (i = 0; jump_string[i] != '\0'; i++) {
 	    bits = jump_string[i];
 	    assert(isxdigit(bits));
 	    bits = tolower(bits);
@@ -164,7 +166,7 @@ extern "C" {
 		bits = bits - '0';
 	    }
 	    bits = bits & 0x0f;
-	    for (int j = 0; j < 4; j++) {
+	    for (j = 0; j < 4; j++) {
 		if ((bits & 1) != 0) {
 		    add(&work, dsfmt);
 		}
diff --git a/setup.py b/setup.py
@@ -124,6 +124,7 @@ def write_config(file_name, config):
 
     elif rng == 'RNG_DSFMT':
         sources += [join(mod_dir, 'src', 'dSFMT', 'dSFMT.c')]
+        sources += [join(mod_dir, 'src', 'dSFMT', 'dSFMT-jump.c')]
         sources += [join(mod_dir, 'shims', 'dSFMT', 'dSFMT-shim.c')]
         # TODO: HAVE_SSE2 should only be for platforms that have SSE2
         # TODO: But how to reliably detect?
