53 * least approximately, for others as well. |
53 * least approximately, for others as well. |
54 * <p> |
54 * <p> |
55 * {@link L128X256MixRandom} is a specific member of the LXM family of algorithms |
55 * {@link L128X256MixRandom} is a specific member of the LXM family of algorithms |
56 * for pseudorandom number generators. Every LXM generator consists of two |
56 * for pseudorandom number generators. Every LXM generator consists of two |
57 * subgenerators; one is an LCG (Linear Congruential Generator) and the other is |
57 * subgenerators; one is an LCG (Linear Congruential Generator) and the other is |
58 * an Xorshift generator. Each output of an LXM generator is the sum of one |
58 * an Xorshift generator. Each output of an LXM generator is the result of |
59 * output from each subgenerator, possibly processed by a final mixing function |
59 * combining state from the LCG with state from the Xorshift generator by |
60 * (and {@link L128X256MixRandom} does use a mixing function). |
60 * using a Mixing function (and then the state of the LCG and the state of the |
|
61 * Xorshift generator are advanced). |
61 * <p> |
62 * <p> |
62 * The LCG subgenerator for {@link L128X256MixRandom} has an update step of the |
63 * The LCG subgenerator for {@link L128X256MixRandom} has an update step of the |
63 * form {@code s = m * s + a}, where {@code s}, {@code m}, and {@code a} are all |
64 * form {@code s = m * s + a}, where {@code s}, {@code m}, and {@code a} are all |
64 * 128-bit integers; {@code s} is the mutable state, the multiplier {@code m} |
65 * 128-bit integers; {@code s} is the mutable state, the multiplier {@code m} |
65 * is fixed (the same for all instances of {@link L128X256MixRandom}) and the addend |
66 * is fixed (the same for all instances of {@link L128X256MixRandom}) and the addend |
72 * version 1.0 (parameters 17, 45), without any final scrambler such as "+" or "**". |
73 * version 1.0 (parameters 17, 45), without any final scrambler such as "+" or "**". |
73 * Its state consists of four {@code long} fields {@code x0}, {@code x1}, {@code x2}, |
74 * Its state consists of four {@code long} fields {@code x0}, {@code x1}, {@code x2}, |
74 * and {@code x3}, which can take on any values provided that they are not all zero. |
75 * and {@code x3}, which can take on any values provided that they are not all zero. |
75 * The period of this subgenerator is 2<sup>256</sup>-1. |
76 * The period of this subgenerator is 2<sup>256</sup>-1. |
76 * <p> |
77 * <p> |
77 * The mixing function for {@link L128X256MixRandom} is the 64-bit MurmurHash3 finalizer. |
78 * The mixing function for {@link L128X256MixRandom} is {@link RandomSupport.mixLea64} |
|
79 * applied to the argument {@code (sh + x0)}, where {@code sh} is the high half of {@code s}. |
78 * <p> |
80 * <p> |
79 * Because the periods 2<sup>128</sup> and 2<sup>256</sup>-1 of the two subgenerators |
81 * Because the periods 2<sup>128</sup> and 2<sup>256</sup>-1 of the two subgenerators |
80 * are relatively prime, the <em>period</em> of any single {@link L128X256MixRandom} object |
82 * are relatively prime, the <em>period</em> of any single {@link L128X256MixRandom} object |
81 * (the length of the series of generated 64-bit values before it repeats) is the product |
83 * (the length of the series of generated 64-bit values before it repeats) is the product |
82 * of the periods of the subgenerators, that is, 2<sup>128</sup>(2<sup>256</sup>-1), |
84 * of the periods of the subgenerators, that is, 2<sup>128</sup>(2<sup>256</sup>-1), |
84 * {@link L128X256MixRandom} objects have different {@code a} parameters, then their |
86 * {@link L128X256MixRandom} objects have different {@code a} parameters, then their |
85 * cycles of produced values will be different. |
87 * cycles of produced values will be different. |
86 * <p> |
88 * <p> |
87 * The 64-bit values produced by the {@code nextLong()} method are exactly equidistributed. |
89 * The 64-bit values produced by the {@code nextLong()} method are exactly equidistributed. |
88 * For any specific instance of {@link L128X256MixRandom}, over the course of its cycle each |
90 * For any specific instance of {@link L128X256MixRandom}, over the course of its cycle each |
89 * of the 2<sup>64</sup> possible {@code long} values will be produced 2<sup>256</sup>-1 times. |
91 * of the 2<sup>64</sup> possible {@code long} values will be produced |
90 * The values produced by the {@code nextInt()}, {@code nextFloat()}, and {@code nextDouble()} |
92 * 2<sup>64</sup>(2<sup>256</sup>-1) times. The values produced by the {@code nextInt()}, |
91 * methods are likewise exactly equidistributed. |
93 * {@code nextFloat()}, and {@code nextDouble()} methods are likewise exactly equidistributed. |
92 * <p> |
94 * <p> |
93 * In fact, the 64-bit values produced by the {@code nextLong()} method are exactly |
95 * Moreover, 64-bit values produced by the {@code nextLong()} method are conjectured to be |
94 * 2-equidistributed. For any specific instance of {@link L128X256MixRandom}, consider |
96 * "very nearly" 4-equidistributed: all possible quadruples of 64-bit values are generated, |
95 * the (overlapping) length-2 subsequences of the cycle of 64-bit values produced by |
97 * and some pairs occur more often than others, but only very slightly more often. |
96 * {@code nextLong()} (assuming no other methods are called that would affect the state). |
98 * However, this conjecture has not yet been proven mathematically. |
97 * There are 2<sup>128</sup>(2<sup>256</sup>-1) such subsequences, and each subsequence, |
99 * If this conjecture is true, then the values produced by the {@code nextInt()}, {@code nextFloat()}, |
98 * which consists of 2 64-bit values, can have one of 2<sup>128</sup> values, and each |
100 * and {@code nextDouble()} methods are likewise approximately 4-equidistributed. |
99 * such value occurs 2<sup>256</sup>-1 times. The values produced by the {@code nextInt()}, |
|
100 * {@code nextFloat()}, and {@code nextDouble()} methods are likewise exactly 2-equidistributed. |
|
101 * <p> |
|
102 * Moreover, the 64-bit values produced by the {@code nextLong()} method are 4-equidistributed. |
|
103 * To be precise: for any specific instance of {@link L128X256MixRandom}, consider |
|
104 * the (overlapping) length-4 subsequences of the cycle of 64-bit values produced by |
|
105 * {@code nextLong()} (assuming no other methods are called that would affect the state). |
|
106 * There are <sup>128</sup>(2<sup>256</sup>-1) such subsequences, and each subsequence, |
|
107 * which consists of 4 64-bit values, can have one of 2<sup>256</sup> values. Of those |
|
108 * 2<sup>256</sup> subsequence values, nearly all of them (2<sup>256</sup>-2<sup>128</sup>) |
|
109 * occur 2<sup>128</sup> times over the course of the entire cycle, and the other |
|
110 * 2<sup>128</sup> subsequence values occur only 2<sup>128</sup>-1 times. So the ratio |
|
111 * of the probability of getting one of the less common subsequence values and the |
|
112 * probability of getting one of the more common subsequence values is 1-2<sup>-128</sup>. |
|
113 * (Note that the set of 2<sup>128</sup> less-common subsequence values will differ from |
|
114 * one instance of {@link L128X256MixRandom} to another, as a function of the additive |
|
115 * parameter of the LCG.) The values produced by the {@code nextInt()}, {@code nextFloat()}, |
|
116 * and {@code nextDouble()} methods are likewise 4-equidistributed. |
|
117 * <p> |
101 * <p> |
118 * Method {@link #split} constructs and returns a new {@link L128X256MixRandom} |
102 * Method {@link #split} constructs and returns a new {@link L128X256MixRandom} |
119 * instance that shares no mutable state with the current instance. However, with |
103 * instance that shares no mutable state with the current instance. However, with |
120 * very high probability, the values collectively generated by the two objects |
104 * very high probability, the values collectively generated by the two objects |
121 * have the same statistical properties as if the same quantity of values were |
105 * have the same statistical properties as if the same quantity of values were |
191 */ |
175 */ |
192 private static final BigInteger PERIOD = |
176 private static final BigInteger PERIOD = |
193 BigInteger.ONE.shiftLeft(256).subtract(BigInteger.ONE).shiftLeft(128); |
177 BigInteger.ONE.shiftLeft(256).subtract(BigInteger.ONE).shiftLeft(128); |
194 |
178 |
195 /* |
179 /* |
196 * The multiplier used in the LCG portion of the algorithm is 2**64 + m; |
180 * Low half of multiplier used in the LCG portion of the algorithm; |
197 * where m is taken from |
181 * the overall multiplier is (2**64 + ML). |
198 * Pierre L'Ecuyer, Tables of linear congruential generators of |
182 * Chosen based on research by Sebastiano Vigna and Guy Steele (2019). |
199 * different sizes and good lattice structure, <em>Mathematics of |
183 * The spectral scores for dimensions 2 through 8 for the multiplier 0x1d605bbb58c8abbfdLL |
200 * Computation</em> 68, 225 (January 1999), pages 249-260, |
184 * are [0.991889, 0.907938, 0.830964, 0.837980, 0.780378, 0.797464, 0.761493]. |
201 * Table 4 (first multiplier for size 2<sup>64</sup>). |
185 */ |
202 * |
186 |
203 * This is almost certainly not the best possible 128-bit multiplier |
187 private static final long ML = 0xd605bbb58c8abbfdL; |
204 * for an LCG, but it is sufficient for our purposes here; because |
|
205 * is is larger than 2**64, the 64-bit values produced by nextLong() |
|
206 * are exactly 2-equidistributed, and the fact that it is of the |
|
207 * form (2**64 + m) simplifies the code, given that we have only |
|
208 * 64-bit arithmetic to work with. |
|
209 */ |
|
210 |
|
211 private static final long M = 2862933555777941757L; |
|
212 |
188 |
213 /* ---------------- instance fields ---------------- */ |
189 /* ---------------- instance fields ---------------- */ |
214 |
190 |
215 /** |
191 /** |
216 * The parameter that is used as an additive constant for the LCG. |
192 * The parameter that is used as an additive constant for the LCG. |
217 * Must be odd. |
193 * Must be odd (therefore al must be odd). |
218 */ |
194 */ |
219 private final long ah, al; |
195 private final long ah, al; |
220 |
196 |
221 /** |
197 /** |
222 * The per-instance state: sh and sl for the LCG; x0, x1, x2, and x3 for the xorshift. |
198 * The per-instance state: sh and sl for the LCG; x0, x1, x2, and x3 for the xorshift. |
250 this.x1 = x1; |
226 this.x1 = x1; |
251 this.x2 = x2; |
227 this.x2 = x2; |
252 this.x3 = x3; |
228 this.x3 = x3; |
253 // If x0, x1, x2, and x3 are all zero, we must choose nonzero values. |
229 // If x0, x1, x2, and x3 are all zero, we must choose nonzero values. |
254 if ((x0 | x1 | x2 | x3) == 0) { |
230 if ((x0 | x1 | x2 | x3) == 0) { |
|
231 long v = sh; |
255 // At least three of the four values generated here will be nonzero. |
232 // At least three of the four values generated here will be nonzero. |
256 this.x0 = RandomSupport.mixStafford13(sh += RandomSupport.GOLDEN_RATIO_64); |
233 this.x0 = RandomSupport.mixStafford13(v += RandomSupport.GOLDEN_RATIO_64); |
257 this.x1 = RandomSupport.mixStafford13(sh += RandomSupport.GOLDEN_RATIO_64); |
234 this.x1 = RandomSupport.mixStafford13(v += RandomSupport.GOLDEN_RATIO_64); |
258 this.x2 = RandomSupport.mixStafford13(sh += RandomSupport.GOLDEN_RATIO_64); |
235 this.x2 = RandomSupport.mixStafford13(v += RandomSupport.GOLDEN_RATIO_64); |
259 this.x3 = RandomSupport.mixStafford13(sh + RandomSupport.GOLDEN_RATIO_64); |
236 this.x3 = RandomSupport.mixStafford13(v + RandomSupport.GOLDEN_RATIO_64); |
260 } |
237 } |
261 } |
238 } |
262 |
239 |
263 /** |
240 /** |
264 * Creates a new instance of {@link L128X256MixRandom} using the |
241 * Creates a new instance of {@link L128X256MixRandom} using the |
275 // |
252 // |
276 // The seed is hashed by mixMurmur64 to produce the `a` parameter. |
253 // The seed is hashed by mixMurmur64 to produce the `a` parameter. |
277 // The seed is hashed by mixStafford13 to produce the initial `x0`, |
254 // The seed is hashed by mixStafford13 to produce the initial `x0`, |
278 // which will then be used to produce the first generated value. |
255 // which will then be used to produce the first generated value. |
279 // The other x values are filled in as if by a SplitMix PRNG with |
256 // The other x values are filled in as if by a SplitMix PRNG with |
280 // GOLDEN_RATIO_64 as the gamma value and Stafford13 as the mixer. |
257 // GOLDEN_RATIO_64 as the gamma value and mixStafford13 as the mixer. |
281 this(RandomSupport.mixMurmur64(seed ^= RandomSupport.SILVER_RATIO_64), |
258 this(RandomSupport.mixMurmur64(seed ^= RandomSupport.SILVER_RATIO_64), |
282 RandomSupport.mixMurmur64(seed += RandomSupport.GOLDEN_RATIO_64), |
259 RandomSupport.mixMurmur64(seed += RandomSupport.GOLDEN_RATIO_64), |
283 0, |
260 0, |
284 1, |
261 1, |
285 RandomSupport.mixStafford13(seed), |
262 RandomSupport.mixStafford13(seed), |
321 this.x2 = x2; |
298 this.x2 = x2; |
322 this.x3 = x3; |
299 this.x3 = x3; |
323 } |
300 } |
324 |
301 |
325 /* ---------------- public methods ---------------- */ |
302 /* ---------------- public methods ---------------- */ |
326 |
303 |
327 /** |
304 /** |
328 * Constructs and returns a new instance of {@link L128X256MixRandom} |
305 * Given 63 bits of "brine", constructs and returns a new instance of |
329 * that shares no mutable state with this instance. |
306 * {@code L128X256MixRandom} that shares no mutable state with this instance. |
330 * However, with very high probability, the set of values collectively |
307 * However, with very high probability, the set of values collectively |
331 * generated by the two objects has the same statistical properties as if |
308 * generated by the two objects has the same statistical properties as if |
332 * same the quantity of values were generated by a single thread using |
309 * same the quantity of values were generated by a single thread using |
333 * a single {@link L128X256MixRandom} object. Either or both of the two |
310 * a single {@code L128X256MixRandom} object. Either or both of the two |
334 * objects may be further split using the {@code split} method, |
311 * objects may be further split using the {@code split} method, |
335 * and the same expected statistical properties apply to the |
312 * and the same expected statistical properties apply to the |
336 * entire set of generators constructed by such recursive splitting. |
313 * entire set of generators constructed by such recursive splitting. |
337 * |
314 * |
338 * @param source a {@link SplittableGenerator} instance to be used instead |
315 * @param source a {@code SplittableGenerator} instance to be used instead |
339 * of this one as a source of pseudorandom bits used to |
316 * of this one as a source of pseudorandom bits used to |
340 * initialize the state of the new ones. |
317 * initialize the state of the new ones. |
341 * @return a new instance of {@link L128X256MixRandom} |
318 * @param brine a long value, of which the low 63 bits are used to choose |
342 */ |
319 * the {@code a} parameter for the new instance. |
343 public L128X256MixRandom split(SplittableGenerator source) { |
320 * @return a new instance of {@code L128X256MixRandom} |
344 // Literally pick a new instance "at random". |
321 */ |
345 return new L128X256MixRandom(source.nextLong(), source.nextLong(), |
322 public SplittableGenerator split(SplittableGenerator source, long brine) { |
346 source.nextLong(), source.nextLong(), |
323 // Pick a new instance "at random", but use the brine for (the low half of) `a`. |
347 source.nextLong(), source.nextLong(), |
324 return new L128X256MixRandom(source.nextLong(), brine << 1, |
348 source.nextLong(), source.nextLong()); |
325 source.nextLong(), source.nextLong(), |
|
326 source.nextLong(), source.nextLong(), |
|
327 source.nextLong(), source.nextLong()); |
349 } |
328 } |
350 |
329 |
351 /** |
330 /** |
352 * Returns a pseudorandom {@code long} value. |
331 * Returns a pseudorandom {@code long} value. |
353 * |
332 * |
354 * @return a pseudorandom {@code long} value |
333 * @return a pseudorandom {@code long} value |
355 */ |
334 */ |
356 public long nextLong() { |
335 public long nextLong() { |
357 final long z = sh + x0; |
336 // Compute the result based on current state information |
358 // The LCG: in effect, s = ((1LL << 64) + M) * s + a, if only we had 128-bit arithmetic. |
337 // (this allows the computation to be overlapped with state update). |
359 final long u = M * sl; |
338 final long result = RandomSupport.mixLea64(sh + x0); |
360 sh = (M * sh) + Math.multiplyHigh(M, sl) + sl + ah; |
339 |
|
340 // Update the LCG subgenerator |
|
341 // The LCG is, in effect, s = ((1LL << 64) + ML) * s + a, if only we had 128-bit arithmetic. |
|
342 final long u = ML * sl; |
|
343 // Note that Math.multiplyHigh computes the high half of the product of signed values, |
|
344 // but what we need is the high half of the product of unsigned values; for this we use the |
|
345 // formula "unsignedMultiplyHigh(a, b) = multiplyHigh(a, b) + ((a >> 63) & b) + ((b >> 63) & a)"; |
|
346 // in effect, each operand is added to the result iff the sign bit of the other operand is 1. |
|
347 // (See Henry S. Warren, Jr., _Hacker's Delight_ (Second Edition), Addison-Wesley (2013), |
|
348 // Section 8-3, p. 175; or see the First Edition, Addison-Wesley (2003), Section 8-3, p. 133.) |
|
349 // If Math.unsignedMultiplyHigh(long, long) is ever implemented, the following line can become: |
|
350 // sh = (ML * sh) + Math.unsignedMultiplyHigh(ML, sl) + sl + ah; |
|
351 // and this entire comment can be deleted. |
|
352 sh = (ML * sh) + (Math.multiplyHigh(ML, sl) + ((ML >> 63) & sl) + ((sl >> 63) & ML)) + sl + ah; |
361 sl = u + al; |
353 sl = u + al; |
362 if (Long.compareUnsigned(sl, u) < 0) ++sh; // Handle the carry propagation from low half to high half. |
354 if (Long.compareUnsigned(sl, u) < 0) ++sh; // Handle the carry propagation from low half to high half. |
|
355 |
|
356 // Update the Xorshift subgenerator |
363 long q0 = x0, q1 = x1, q2 = x2, q3 = x3; |
357 long q0 = x0, q1 = x1, q2 = x2, q3 = x3; |
364 { // xoshiro256 1.0 |
358 { // xoshiro256 1.0 |
365 long t = q1 << 17; |
359 long t = q1 << 17; |
366 q2 ^= q0; |
360 q2 ^= q0; |
367 q3 ^= q1; |
361 q3 ^= q1; |