Skip to content

Commit 4d4228e

Browse files
bashtagebashtage
authored
Merge pull request #105 from bashtage/sfmt-jump
ENH: Add jump for sfmt
2 parents 2f66258 + 713fa6e commit 4d4228e

File tree

15 files changed

+690
-26
lines changed

15 files changed

+690
-26
lines changed

‎doc/source/change-log.rst

Lines changed: 8 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -3,8 +3,8 @@
33
Change Log
44
==========
55

6-
Since Version 1.13
7-
------------------
6+
Version 1.13.2
7+
--------------
88
* Add Ziggurat generation for standard gamma
99
(:meth:`~randomstate.prng.mt19937.standard_gamma`) for both floats and
1010
doubles. The gamma generator uses a rejection sampler that
@@ -34,6 +34,7 @@ Since Version 1.13
3434
(`SFMT <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/>`_) generator.
3535
* Add complex normal (:meth:`~randomstate.prng.mt19937.complex_normal`)
3636
* Added support for jumping the MT19937 generator
37+
* Added support for jumping the SFMT generator
3738

3839
Version 1.13
3940
------------
@@ -59,8 +60,10 @@ Version 1.11.3
5960
--------------
6061
* Extended 32-bit generation to
6162

62-
* Uniforms (:meth:`~randomstate.prng.mt19937.random_sample` and :meth:`~randomstate.prng.mt19937.rand`)
63-
* Normals (:meth:`~randomstate.prng.mt19937.standard_normal` and :meth:`~randomstate.prng.mt19937.randn`)
63+
* Uniforms (:meth:`~randomstate.prng.mt19937.random_sample` and
64+
:meth:`~randomstate.prng.mt19937.rand`)
65+
* Normals (:meth:`~randomstate.prng.mt19937.standard_normal` and
66+
:meth:`~randomstate.prng.mt19937.randn`)
6467
* Standard Gammas (:meth:`~randomstate.prng.mt19937.standard_gamma`)
6568
* Standard Exponentials (:meth:`~randomstate.prng.mt19937.standard_exponential`)
6669

@@ -95,7 +98,7 @@ Version 1.11.1
9598
Version 1.11
9699
------------
97100
* Update to recent changes in NumPy's RandomState
98-
* Expose system entropy through :meth:`randomstate.entropy.random_entropy`
101+
* Expose system entropy through :func:`randomstate.entropy.random_entropy`
99102
* Add vector initialization for all PRNGs
100103

101104
Version 1.10.1

‎doc/source/entropy.rst

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,6 @@
11
System Entropy
22
==============
33

4-
.. autofunction:: randomstate.entropy.random_entropy
4+
.. module:: randomstate.entropy
5+
6+
.. autofunction:: random_entropy

‎doc/source/index.rst

Lines changed: 18 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -5,14 +5,16 @@ that change the core random number generator.
55

66
What's New or Different
77
-----------------------
8-
* ``randomstate.entropy.random_entropy`` provides access to the system
8+
* :func:`randomstate.entropy.random_entropy` provides access to the system
99
source of randomness that is used in cryptographic applications (e.g.,
1010
``/dev/urandom`` on Unix).
1111
* Simulate from the complex normal distribution
1212
(:meth:`~randomstate.prng.mt19937.complex_normal`)
1313
* The normal, exponential and gamma generators support 256-step Ziggurat
1414
methods which are 2-10 times faster than NumPy's default implementation in
15-
``standard_normal``, ``standard_exponential`` or ``standard_gamma``.
15+
:meth:`~randomstate.prng.mt19937.standard_normal`,
16+
:meth:`~randomstate.prng.mt19937.standard_exponential` or
17+
:meth:`~randomstate.prng.mt19937.standard_gamma`.
1618
The Ziggurat generator can be accessed by passing the keyword
1719
argument ``method='zig'``.
1820

@@ -38,8 +40,10 @@ What's New or Different
3840
to produce either single or double prevision uniform random variables for
3941
select distributions
4042

41-
* Uniforms (:meth:`~randomstate.prng.mt19937.random_sample` and :meth:`~randomstate.prng.mt19937.rand`)
42-
* Normals (:meth:`~randomstate.prng.mt19937.standard_normal` and :meth:`~randomstate.prng.mt19937.randn`)
43+
* Uniforms (:meth:`~randomstate.prng.mt19937.random_sample` and
44+
:meth:`~randomstate.prng.mt19937.rand`)
45+
* Normals (:meth:`~randomstate.prng.mt19937.standard_normal` and
46+
:meth:`~randomstate.prng.mt19937.randn`)
4347
* Standard Gammas (:meth:`~randomstate.prng.mt19937.standard_gamma`)
4448
* Standard Exponentials (:meth:`~randomstate.prng.mt19937.standard_exponential`)
4549

@@ -91,16 +95,18 @@ generators, 'in addition' to the standard PRNG in NumPy. The included PRNGs are
9195
using the same seed/state. Adds a jump function that advances the generator
9296
as-if 2**128 draws have been made (:meth:`randomstate.prng.mt19937.jump`).
9397
See `NumPy's documentation`_.
94-
* dSFMT - A SSE2 enables version of the MT19937 generator. Theoretically the
95-
same, but with a different state and so it is not possible to produce a
96-
sequence identical to MT19937. See the `dSFMT authors' page`_.
97-
* XoroShiro128+ - Improved version of XorShift128+ with improved performance
98-
and better statistical quality. Like the XorShift generators, can be jumped
98+
* SFMT and dSFMT - SSE2 enabled versions of the MT19937 generator. Theoretically
99+
the same, but with a different state and so it is not possible to produce a
100+
sequence identical to MT19937. Both generators support ``jump`` and so can
101+
be used in parallel applications. See the `dSFMT authors' page`_.
102+
* XoroShiro128+ - Improved version of XorShift128+ with better performance
103+
and statistical quality. Like the XorShift generators, it can be jumped
99104
to produce multiple streams in parallel applications. See
100-
:meth:`randomstate.prng.xoroshiro128plus.jump` for details. More information
101-
about this PRNG is available at the `xorshift and xoroshiro authors' page`_.
105+
:meth:`randomstate.prng.xoroshiro128plus.jump` for details.
106+
More information about this PRNG is available at the
107+
`xorshift and xoroshiro authors' page`_.
102108
* XorShit128+ and XorShift1024* - Vast fast generators based on the XSadd
103-
generator. These generators can be rapidly 'jumped' and so can be used in
109+
generator. These generators support ``jump`` and so can be used in
104110
parallel applications. See the documentation for
105111
:meth:`randomstate.prng.xorshift1024.jump` for details. More information
106112
about these PRNGs is available at the

‎doc/source/sfmt.rst

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,13 @@ Random generator
1515
get_state
1616
set_state
1717

18+
Parallel generation
19+
===================
20+
.. autosummary::
21+
:toctree: generated/
22+
23+
jump
24+
1825
Simple random data
1926
==================
2027
.. autosummary::

‎randomstate/interface/sfmt/sfmt-poly.h

Lines changed: 69 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,69 @@
1+
static const char *poly_128 =
2+
"c5ecf0b605bcbebf05a6ae430dc3cd1777acaff2b99fe5221eaa2a2c8acf0ea0b9617d0b4a"
3+
"5e8f94828b6710c0f365e69b690be81b206e04f64a9994b013dc1e046fb5e182b3d794fdd3"
4+
"d01d73d1c2aaaa03ce4398cce92cc67ed20c031a607fdca6dce456b8268481c80bcf740a81"
5+
"0c5d73dcab943e05e8a2adfe5f81d9f118435ceef77fbe3d9d58baa17378e24001fea40b87"
6+
"181de63568f4342c05ddd174a617d1452cfcb4fd1df0dc94802ce7b8c04f21ced8cd34d297"
7+
"71f5cb25471704d2bfc8abc2acc47f10d9f72f781556c62e4a089f44c63d4a2d7b3541fe07"
8+
"c197c0bd10178c2a81fa7888d0ddc777502f8dda122e3bf335e365f0e70be77f0383bf2d0b"
9+
"5ba5c1ab9cd94288e59f086262c42bde8f832819e9fa1467661d68ff8ed3edd81c8aa6c268"
10+
"a72bcafbc3ab20c1ecb1799988abf25a6b9b4262bf93fbd2ef44fe16e6a2702e57b8d1f046"
11+
"8a58a85cc74f2a75d18b3373cfc47cfd028f70ea8e6a5ef2e871b929b47f4fd1ffe3be9c4c"
12+
"89ac6c2c135fd1f410c8e0e9e55e4a75891f00e3a48f6a4a48b1e12f9032675e4942c50224"
13+
"1e3d75b34fe803d6cd4e33709cb2632a4fe42c7a94d1a57b166f1a9d19e58bc7d904740df1"
14+
"bca2481a01e6f390d647353dd18dc1ba77ac5334e3e0037413f2f41d33941ecfdd1cbdf0cb"
15+
"e4ad3c08ac47ec03c440971327b1631a8742ef3fa2d2fa7e76af3de0216f52d7f44d0efcaf"
16+
"927f54b673a0ade6910f7fb86da9ef4f3715d84fcd2159b3133371bdbf6e888276e238e3b0"
17+
"3ee0c255000ae8914736f597f7222c1764eb1b13a9447c406257b23e707576b0456b76fc41"
18+
"4a75620093e8f7017ac980a867f1221c1eb0979d14032c0838d976cf5d139b1c5a610ab09f"
19+
"117776c2167b201056b7f0cc246232d183c207ae38dd8d0fccc302b034bcdeb4976b2b8864"
20+
"cc5fc2f162b72b55f3ee9cd28fd4f271063d4a58e3a66698d56487651dbf7a135b383adb71"
21+
"661d188f6c39521381dbc3b2bd01041f2ebd396f9f76511cffc7ef7897903bdee6c086e193"
22+
"c541ebb401dca2d28d779000a2d69740e8f6560477c23b01920175c38604150200ecce4d4c"
23+
"cfd0a6cb6f4f2bf2ab37a3072d4e613356307338ee2e6767b6a83ec6cda267bb1994e099af"
24+
"4067d67d95341a1ca150f7196db5e2209969b3d432179a3d8da93f71825bc7b385d960f4b8"
25+
"81b30aa6a4cf16a2672eda09b3efea2ed7cbcf77ed405f1aa175c59f5f619a54cc0a437ba9"
26+
"7cd8d041a5ca0e75534b2e9c61d5ee2e1570945fe276921e18308db9f3a9805d798611f7f4"
27+
"96d2b958b6eb4e946d572c471c619565832d20d0e82eae1eef9e1a2472c6c851250c08aae4"
28+
"d403558e1fde9c399e324a5e8dede8ca57b4492afda60cea0996508a3567a01a5e6e822bca"
29+
"9ec8ac223bc0534c311ed07eae4bcaed013776ba0bbdb22b8bc038b61cb431c8079201e0f8"
30+
"304a439fb8dc1d71050e4d6c93b154a152def253ad808a11dea33c6d797c94ac0eac9b5180"
31+
"96ccb85081c2c02a89e32624e700cd1ccfa7eb90592afaf886d4f61f1a378ab44f29dd73df"
32+
"e5ab7a23a8f31a9b6f12639c81d3db851abff93938a60f61b659a481f286ba10339e0409f3"
33+
"9885f58b42e255ad89d0484763459e03a04097d2427444f59e3d49206f784ef676914f4a4e"
34+
"11b66635881ab05accd0020559a780ebd53999974dae64cfedc4ee856ae2d3b3cd5d40337d"
35+
"33115dea8376f07da698f20c11dc86d8203d0d12bf3ae973bfec0700d2337961e1c4c9448c"
36+
"b5730d9457afb1d8891301f3264fd30cba095db3f527907c80fca305ea8b9ca2adcdb0e2b2"
37+
"781839a1ebcb7ff5b627ddfcf72170118606e5742be995af9d5e0f31798640f32d6e09438c"
38+
"c75737d6876daa5c0eea191c6aceff45d66b5beceda17833b2968b1e9454612f85a93c239c"
39+
"96a96eead57a36da005c661e1ac83427e5724b06b7423e99731bd126719998e17ab3be9d2c"
40+
"141002533c84cd27ff029cea709f5c585dc1194bcb6df9663a7b1c6adca72895bc2f42bf80"
41+
"19f52975398d8d68e64b76796e93187ebb1d92a10387c5d41687d8bacebf2705904783edd7"
42+
"1b1b652f4acf1185289441b21ddb079de827dad3adfc68a1832a2ef50d9b4b324392267dd3"
43+
"ed9094eb442c4306758ba07225268026010bf582ce2563dc4de2b6ad0c8b0aac39003f4b0e"
44+
"c473854125038a77be70ec8d41301c56e970d8aef57ce13f0619f882edc582d78248178adf"
45+
"5576466b0215ae8b899d95df9b6753c87fdc3a54ed8a90923dce88b89f371235ca1d0f4c0d"
46+
"6737ee8db168f7ef4aea9c3bdf06d9e96236333b7ee207bdf4c3603c55be04a890837a041f"
47+
"3a8e6ec6e3f91e57a94504098e341ec44a3236608d439cd0b5840d7f1f2c589fb356f378c6"
48+
"6439dc97ee5b176050a208bc813a2db4b023dbc5298a013aefc56792d62c245ff9d18228d5"
49+
"3ae155918d0a6ccf61131cc5cdf34deb4545eb48a0f5b7e486b146af6699a756fac06e9a33"
50+
"1a8fa944a5b0046321a40635a709bf2b9bbff44d907d83924dc0d908ee1d21c60100ddb448"
51+
"687928ed760783568ec89d65fdbd077d5312683ca526d413f4bee45e814a160854b7329c83"
52+
"fe490a242595fbf3b85f32b4715f1034c9737d7fab25b4b29b4c90c9919df8377c5111cfdb"
53+
"b43e6ad9e467ff48d0d85de6457131a4b4ca6a69e3de7d09691a9b5da24dda54f3fdc5f5c3"
54+
"16ddcc99842cce8af4440034f4bf6236fae7092d41a027880b8e29371a476451997d52ec9b"
55+
"c428308e171499e29fd008f7a9dceff044deb0e7164e774b248e472ef5f38893e7f71bfbfb"
56+
"4728e6063ba1111df7b978e81313388f5b6bc1111ca35885810ed57cc6fce9a8f53c7cf40c"
57+
"c360ef836bc73a937a1979117466b2e1ab4a11cbc0be9c0741d12d5d1f7c7333f7d34ec9d8"
58+
"734541f0544a7242895a40eec2d46e606d751b275ac6b7d208bb5cd112dedea4652a38b134"
59+
"ad061771b9ba842108e77ad0149d3fb265a277a770d9083825b7876e8c30a1f83dec6d54da"
60+
"4da61dbeadf706e6a8d18d4e00fcd19985cb61c3e25439cfa68ba1f007b43a4d79c1f2ae85"
61+
"66318e0d0c0ac6e312fd9af2760258a16c1c2e864ec053ce4970520e1780c9ffa218fcbeeb"
62+
"803e553f29ed66da093b8498ca20e030c9dab6d2ab214837468cc65c3a122222d482917c2e"
63+
"f4aae56e8b0f7625997cc5e726ae9b110ab111b123e71e123ba0b2847481d324d6f5568878"
64+
"613ab22aa7fe23ac40b04c82d9c118803804aa014d7b60d1d05c5bd4ff07e44cdedf09f8db"
65+
"d3738366d87301460db172acba11fcdba924ed7bc0debeabc4ea9ed658d12abc8d529f8b9b"
66+
"c736b07968db946a80c9a6e0d5ba9b412742aaa0675b09844685059cd00ff7ab088813ed6d"
67+
"b4460004d48f61f905c373d223a7363bcc5bbcb6d02bc711ac954b0a694c7bc40dbbb9eeef"
68+
"50282e63c56088786a70b1fc8720aeb993451c33a533e3703ce76e2eb8450796109cff8c2d"
69+
"4ba3f2e1dcd6d9397608ac967baa41b564";

‎randomstate/interface/sfmt/sfmt-shim.c

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
11
#include "sfmt-shim.h"
2+
#include "sfmt-poly.h"
23

34
extern inline uint32_t random_uint32(aug_state *state);
45

@@ -33,7 +34,7 @@ void entropy_init(aug_state *state) {
3334
set_seed(state, seeds[0]);
3435
}
3536

36-
// void jump_state(aug_state* state)
37-
//{
38-
// dSFMT_jump(state->rng, poly_128);
39-
//}
37+
void jump_state(aug_state* state)
38+
{
39+
SFMT_jump(state->rng, poly_128);
40+
}

‎randomstate/interface/sfmt/sfmt-shim.h

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -65,4 +65,4 @@ extern void set_seed_by_array(aug_state *state, uint32_t init_key[],
6565

6666
extern void set_seed(aug_state *state, uint32_t seed);
6767

68-
// extern void jump_state(aug_state* state);
68+
extern void jump_state(aug_state* state);

‎randomstate/interface/sfmt/sfmt.pxi

Lines changed: 40 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
11
DEF RS_RNG_NAME = u'sfmt'
2+
DEF RS_RNG_JUMPABLE = 1
23
DEF SFMT_MEXP = 19937
34
DEF SFMT_N = 156 # (SFMT_MEXP / 128 + 1)
45
DEF RS_SEED_NBYTES = 1
@@ -121,6 +122,23 @@ number of probability distributions to choose from.
121122
>>> seed = random_entropy()
122123
>>> rs = [rnd.RandomState(seed) for _ in range(10)]
123124
125+
**Parallel Features**
126+
127+
``sfmt.RandomState`` can be used in parallel applications by
128+
calling the method ``jump`` which advances the state as-if :math:`2^{128}`
129+
random numbers have been generated [2]_. This allows the original sequence to
130+
be split so that distinct segments can be used in each worker process. All
131+
generators should be initialized with the same seed to ensure that the
132+
segments come from the same sequence.
133+
134+
>>> from randomstate.entropy import random_entropy
135+
>>> import randomstate.prng.sfmt as rnd
136+
>>> seed = random_entropy()
137+
>>> rs = [rnd.RandomState(seed) for _ in range(10)]
138+
# Advance rs[i] by i jumps
139+
>>> for i in range(10):
140+
rs[i].jump(i)
141+
124142
**State and Seeding**
125143
126144
The ``sfmt.RandomState`` state vector consists of a 624 element array of
@@ -140,3 +158,25 @@ state values.
140158
Twister: a 128-bit Pseudorandom Number Generator." Monte Carlo
141159
and Quasi-Monte Carlo Methods 2006, Springer, pp. 607 -- 622, 2008.
142160
"""
161+
162+
DEF JUMP_DOCSTRING = u"""
163+
jump(iter = 1)
164+
165+
Jumps the state of the random number generator as-if 2**128 random numbers
166+
have been generated.
167+
168+
Parameters
169+
----------
170+
iter : integer, positive
171+
Number of times to jump the state of the prng.
172+
173+
Returns
174+
-------
175+
out : None
176+
Returns 'None' on success.
177+
178+
Notes
179+
-----
180+
Jumping the rng state resets any pre-computed random numbers. This is required
181+
to ensure exact reproducibility.
182+
"""

‎randomstate/src/dSFMT/calc-jump.cpp

Lines changed: 81 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,81 @@
1+
/**
2+
* @file calc-jump.cpp
3+
*
4+
* @brief calc jump function.
5+
*
6+
* @author Mutsuo Saito (Hiroshima University)
7+
* @author Makoto Matsumoto (The University of Tokyo)
8+
*
9+
* Copyright (C) 2012 Mutsuo Saito, Makoto Matsumoto,
10+
* Hiroshima University and The University of Tokyo.
11+
* All rights reserved.
12+
*
13+
* The 3-clause BSD License is applied to this software, see
14+
* LICENSE.txt
15+
*
16+
* Compile:
17+
* g++ calc-jump.cpp -o calc-jump -lntl
18+
*
19+
* Compute polynomial for 2^128 steps:
20+
* ./calc-jump 340282366920938463463374607431768211456 poly.19937.txt
21+
*
22+
*/
23+
#include <iostream>
24+
#include <fstream>
25+
#include <iomanip>
26+
#include <sstream>
27+
#include <string>
28+
#include <inttypes.h>
29+
#include <stdint.h>
30+
#include <time.h>
31+
#include <NTL/GF2X.h>
32+
#include <NTL/vec_GF2.h>
33+
#include <NTL/ZZ.h>
34+
#include "dSFMT-calc-jump.hpp"
35+
36+
using namespace NTL;
37+
using namespace std;
38+
using namespace dsfmt;
39+
40+
static void read_file(GF2X& lcmpoly, long line_no, const string& file);
41+
42+
int main(int argc, char * argv[]) {
43+
if (argc <= 2) {
44+
cout << argv[0] << " jump-step poly-file" << endl;
45+
cout << " jump-step: a number between zero and 2^{DSFMT_MEXP}-1.\n"
46+
<< " large decimal number is allowed." << endl;
47+
cout << " poly-file: one of poly.{MEXP}.txt "
48+
<< "file" << endl;
49+
return -1;
50+
}
51+
string step_string = argv[1];
52+
string filename = argv[2];
53+
long no = 0;
54+
GF2X lcmpoly;
55+
read_file(lcmpoly, no, filename);
56+
ZZ step;
57+
stringstream ss(step_string);
58+
ss >> step;
59+
string jump_str;
60+
calc_jump(jump_str, step, lcmpoly);
61+
cout << "jump polynomial:" << endl;
62+
cout << jump_str << endl;
63+
return 0;
64+
}
65+
66+
67+
static void read_file(GF2X& lcmpoly, long line_no, const string& file)
68+
{
69+
ifstream ifs(file.c_str());
70+
string line;
71+
for (int i = 0; i < line_no; i++) {
72+
ifs >> line;
73+
ifs >> line;
74+
}
75+
if (ifs) {
76+
ifs >> line;
77+
line = "";
78+
ifs >> line;
79+
}
80+
stringtopoly(lcmpoly, line);
81+
}

0 commit comments

Comments
 (0)