-
Notifications
You must be signed in to change notification settings - Fork 2
/
MRG32k3a.java
167 lines (154 loc) · 6.76 KB
/
MRG32k3a.java
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
/* Written in 2019 by Sebastiano Vigna (vigna@acm.org)
To the extent possible under law, the author has dedicated all copyright
and related and neighboring rights to this software to the public domain
worldwide. This software is distributed without any warranty.
See <http://creativecommons.org/publicdomain/zero/1.0/>. */
/**
* A fast, testless implementation based on 64-bit integers of Pierre
* L'Ecuyer's pseudorandom number generator <a href=
* "https://pubsonline.informs.org/doi/abs/10.1287/opre.47.1.159">MRG32k3a</a>.
*
* <p>
* There are three tests in the standard implementation of MRG32k3a: one test to
* correct the combined output, and two tests to correct negative modular
* residuals. In this implementation, the first test is avoided by
* arithmetization, and the other two by ensuring that the argument to the
* modulo operator is nonnegative, which is possible because of the small size
* of the numbers involved.
*
* <p>
* Additionally, the output value is computed on the <em>current</em> state,
* rather then on the next state, to let the processor parallelize internally
* the computation of the output value and of the next state.
*
* <p>
* On an Intel® Core™ i7-7700 CPU @3.60GHz, using Java 13 and
* <a href="http://openjdk.java.net/projects/code-tools/jmh/">JMH</a>
* microbenchmarks the <a href=
* "https://github.com/umontreal-simul/ssj/blob/master/src/main/java/umontreal/ssj/rng/MRG32k3a.java">double-based
* implementation</a> takes ≈21.6 ns to emit a double, the <a href=
* "https://github.com/umontreal-simul/ssj/blob/master/src/main/java/umontreal/ssj/rng/MRG32k3aL.java">official
* long-based implementation</a> takes ≈13.6 ns whereas this implementation
* needs just ≈5.4 ns. (The JMH timings were decreased by 1ns, as using
* the low-level {@code perfasm} profiler the JMH overhead was estimated at
* ≈1ns per call.)
*
* <p>
* The stream generated by this implementation is identical to that of the
* original one, provided that the initial state is set using
* {@linkplain #setSeed(long[]) the provided method} (or the equivalent
* {@linkplain #MRG32k3a(long[]) constructor}). An {@linkplain #setSeed(long)
* additional method} (and {@linkplain #MRG32k3a(long) constructor}) provide a
* simpler way to set the initial state starting from a 64-bit seed.
*
* <p>
* Note that currently this implementation does not offer multiple streams, but
* they could be implemented easily (e.g., borrowing code from the <a href=
* "https://github.com/umontreal-simul/ssj/blob/master/src/main/java/umontreal/ssj/rng/MRG32k3a.java">official
* implementation</a>).
*/
public class MRG32k3a {
private static final long m1 = 4294967087L;
private static final long m2 = 4294944443L;
private static final long a12 = 1403580;
private static final long a13 = 810728;
private static final long a21 = 527612;
private static final long a23 = 1370589;
private static final long corr1 = (m1 * a13);
private static final long corr2 = (m2 * a23);
// This is equivalent to 2.328306549295727688e-10, but more precise as
// there's no rounding involved.
private static final double norm = 0x1.000000d00000bp-32;
private long s10, s11, s12, s20, s21, s22;
private static long staffordMix13(long z) {
z = (z ^ (z >>> 30)) * 0xBF58476D1CE4E5B9L;
z = (z ^ (z >>> 27)) * 0x94D049BB133111EBL;
// This guarantees that we return a positive number
return (z >>> 1) ^ (z >>> 32);
}
/**
* Creates a new generator starting from a 64-bit seed.
*
* @param seed
* a 64-bit integer that will be used to seed a SplitMix64
* generator, whose output will be used to seed this generator.
*/
MRG32k3a(final long seed) {
setSeed(seed);
}
/**
* Creates a new generator starting from a specification of its internal
* state.
*
* @param seed
* six numbers: the first three values (s10, s11, s12; not all
* zeros) must be in the range [0..4294967087) and the last three
* values (s20, s21, s22; not all zeros) in the range
* [0..4294944443).
*/
MRG32k3a(final long... seed) {
setSeed(seed);
}
/**
* Sets the state of this generator starting from a 64-bit seed.
*
* @param seed
* a 64-bit integer that will be used to seed a SplitMix64
* generator, whose output will be used to seed this generator.
*/
public void setSeed(long seed) {
s10 = staffordMix13(seed += 0x9e3779b97f4a7c15L) % m1;
s11 = staffordMix13(seed += 0x9e3779b97f4a7c15L) % m1;
s12 = staffordMix13(seed += 0x9e3779b97f4a7c15L) % m1;
s20 = staffordMix13(seed += 0x9e3779b97f4a7c15L) % m2;
s21 = staffordMix13(seed += 0x9e3779b97f4a7c15L) % m2;
s22 = staffordMix13(seed += 0x9e3779b97f4a7c15L) % m2;
}
/**
* Set the state of this generator.
*
* @param seed
* six numbers: the first three values (s10, s11, s12; not all
* zeros) must be in the range [0..4294967087) and the last three
* values (s20, s21, s22; not all zeros) in the range
* [0..4294944443).
*/
public void setSeed(final long... seed) {
if (seed.length != 6) throw new IllegalArgumentException("You must provide six numbers");
if (seed[0] == 0 && seed[1] == 0 && seed[2] == 0) throw new IllegalArgumentException("s10, s11 and s12 cannot be all zero");
if (seed[3] == 0 && seed[4] == 0 && seed[5] == 0) throw new IllegalArgumentException("s20, s21 and s22 cannot be all zero");
if (seed[0] < 0 || seed[1] < 0 || seed[2] < 0) throw new IllegalArgumentException("s10 (" + seed[0] + "), s11 (" + seed[1] + "), and s12 (" + seed[2] + ") cannot be negative");
if (seed[3] < 0 || seed[4] < 0 || seed[5] < 0) throw new IllegalArgumentException("s20 (" + seed[3] + "), s21 (" + seed[4] + "), and s22 (" + seed[5] + ") cannot be negative");
if (seed[0] >= m1 || seed[1] >= m1 || seed[2] >= m1) throw new IllegalArgumentException("s10 (" + seed[0] + "), s11 (" + seed[1] + "), and s12 (" + seed[2] + ") must be smaller than " + m1);
if (seed[3] >= m2 || seed[4] >= m2 || seed[5] >= m2) throw new IllegalArgumentException("s20 (" + seed[3] + "), s21 (" + seed[4] + "), and s22 (" + seed[5] + ") must be smaller than " + m2);
s10 = seed[0];
s11 = seed[1];
s12 = seed[2];
s20 = seed[3];
s21 = seed[4];
s22 = seed[5];
// Throw away one value to align output with L'Ecuyer original version
nextDouble();
}
/**
* Returns the next pseudorandom double in (0..1).
*
* @return the next pseudorandom double in (0..1).
*/
public double nextDouble() {
// Combination
long r = s12 - s22;
r -= m1 * ((r - 1) >> 63);
// Component 1
long p = (a12 * s11 - a13 * s10 + corr1) % m1;
s10 = s11;
s11 = s12;
s12 = p;
// Component 2
p = (a21 * s22 - a23 * s20 + corr2) % m2;
s20 = s21;
s21 = s22;
s22 = p;
return r * norm;
}
}