20 * |
20 * |
21 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA |
21 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA |
22 * or visit www.oracle.com if you need additional information or have any |
22 * or visit www.oracle.com if you need additional information or have any |
23 * questions. |
23 * questions. |
24 */ |
24 */ |
25 package java.util; |
25 |
26 |
26 package java.util.random; |
27 import java.util.Rng; |
27 |
28 import java.util.Spliterator; |
28 import java.util.Spliterator; |
29 import java.util.function.Consumer; |
29 import java.util.stream.DoubleStream; |
30 import java.util.function.IntConsumer; |
|
31 import java.util.function.LongConsumer; |
|
32 import java.util.function.DoubleConsumer; |
|
33 import java.util.stream.StreamSupport; |
|
34 import java.util.stream.IntStream; |
30 import java.util.stream.IntStream; |
35 import java.util.stream.LongStream; |
31 import java.util.stream.LongStream; |
36 import java.util.stream.DoubleStream; |
|
37 import java.util.DoubleZigguratTables; |
|
38 |
32 |
39 /** |
33 /** |
40 * Low-level utility methods helpful for implementing pseudorandom number generators. |
34 * Low-level utility methods helpful for implementing pseudorandom number generators. |
41 * |
35 * |
42 * This class is mostly for library writers creating specific implementations of the interface {@link java.util.Rng}. |
36 * This class is mostly for library writers creating specific implementations of the |
|
37 * interface {@link RandomNumberGenerator}. |
43 * |
38 * |
44 * @author Guy Steele |
39 * @since 14 |
45 * @author Doug Lea |
|
46 * @since 1.9 |
|
47 */ |
40 */ |
48 public class RngSupport { |
41 public class RNGSupport { |
49 |
42 |
50 /* |
43 /* |
51 * Implementation Overview. |
44 * Implementation Overview. |
52 * |
45 * |
53 * This class provides utility methods and constants frequently |
46 * This class provides utility methods and constants frequently |
54 * useful in the implentation of pseudorandom number generators |
47 * useful in the implentation of pseudorandom number generators |
55 * that satisfy the interface {@code java.util.Rng}. |
48 * that satisfy the interface {@link RandomNumberGenerator}. |
56 * |
49 * |
57 * File organization: First some message strings, then the main |
50 * File organization: First some message strings, then the main |
58 * public methods, followed by a non-public base spliterator class. |
51 * public methods, followed by a non-public base spliterator class. |
59 */ |
52 */ |
60 |
53 |
61 // IllegalArgumentException messages |
54 // IllegalArgumentException messages |
62 static final String BadSize = "size must be non-negative"; |
55 static final String BAD_SIZE = "size must be non-negative"; |
63 static final String BadDistance = "jump distance must be finite, positive, and an exact integer"; |
56 static final String BAD_DISTANCE = "jump distance must be finite, positive, and an exact integer"; |
64 static final String BadBound = "bound must be positive"; |
57 static final String BAD_BOUND = "bound must be positive"; |
65 static final String BadFloatingBound = "bound must be finite and positive"; |
58 static final String BAD_FLOATING_BOUND = "bound must be finite and positive"; |
66 static final String BadRange = "bound must be greater than origin"; |
59 static final String BAD_RANGE = "bound must be greater than origin"; |
67 |
60 |
68 /* ---------------- public methods ---------------- */ |
61 /* ---------------- public methods ---------------- */ |
69 |
62 |
70 /** |
63 /** |
71 * Check a {@code long} proposed stream size for validity. |
64 * Check a {@code long} proposed stream size for validity. |
72 * |
65 * |
73 * @param streamSize the proposed stream size |
66 * @param streamSize the proposed stream size |
|
67 * |
74 * @throws IllegalArgumentException if {@code streamSize} is negative |
68 * @throws IllegalArgumentException if {@code streamSize} is negative |
75 */ |
69 */ |
76 public static void checkStreamSize(long streamSize) { |
70 public static void checkStreamSize(long streamSize) { |
77 if (streamSize < 0L) |
71 if (streamSize < 0L) |
78 throw new IllegalArgumentException(BadSize); |
72 throw new IllegalArgumentException(BAD_SIZE); |
79 } |
73 } |
80 |
74 |
81 /** |
75 /** |
82 * Check a {@code double} proposed jump distance for validity. |
76 * Check a {@code double} proposed jump distance for validity. |
83 * |
77 * |
84 * @param distance the proposed jump distance |
78 * @param distance the proposed jump distance |
85 * @throws IllegalArgumentException if {@code size} not positive, |
79 * |
86 * finite, and an exact integer |
80 * @throws IllegalArgumentException if {@code size} not positive, finite, and an exact integer |
87 */ |
81 */ |
88 public static void checkJumpDistance(double distance) { |
82 public static void checkJumpDistance(double distance) { |
89 if (!(distance > 0.0 && distance < Float.POSITIVE_INFINITY && distance == Math.floor(distance))) |
83 if (!(distance > 0.0 && distance < Float.POSITIVE_INFINITY |
90 throw new IllegalArgumentException(BadDistance); |
84 && distance == Math.floor(distance))) { |
|
85 throw new IllegalArgumentException(BAD_DISTANCE); |
|
86 } |
91 } |
87 } |
92 |
88 |
93 /** |
89 /** |
94 * Checks a {@code float} upper bound value for validity. |
90 * Checks a {@code float} upper bound value for validity. |
95 * |
91 * |
96 * @param bound the upper bound (exclusive) |
92 * @param bound the upper bound (exclusive) |
97 * @throws IllegalArgumentException if {@code bound} is not |
93 * |
98 * positive and finite |
94 * @throws IllegalArgumentException if {@code bound} is not positive and finite |
99 */ |
95 */ |
100 public static void checkBound(float bound) { |
96 public static void checkBound(float bound) { |
101 if (!(bound > 0.0 && bound < Float.POSITIVE_INFINITY)) |
97 if (!(bound > 0.0 && bound < Float.POSITIVE_INFINITY)) { |
102 throw new IllegalArgumentException(BadFloatingBound); |
98 throw new IllegalArgumentException(BAD_FLOATING_BOUND); |
|
99 } |
103 } |
100 } |
104 |
101 |
105 /** |
102 /** |
106 * Checks a {@code double} upper bound value for validity. |
103 * Checks a {@code double} upper bound value for validity. |
107 * |
104 * |
108 * @param bound the upper bound (exclusive) |
105 * @param bound the upper bound (exclusive) |
109 * @throws IllegalArgumentException if {@code bound} is not |
106 * |
110 * positive and finite |
107 * @throws IllegalArgumentException if {@code bound} is not positive and finite |
111 */ |
108 */ |
112 public static void checkBound(double bound) { |
109 public static void checkBound(double bound) { |
113 if (!(bound > 0.0 && bound < Double.POSITIVE_INFINITY)) |
110 if (!(bound > 0.0 && bound < Double.POSITIVE_INFINITY)) { |
114 throw new IllegalArgumentException(BadFloatingBound); |
111 throw new IllegalArgumentException(BAD_FLOATING_BOUND); |
|
112 } |
115 } |
113 } |
116 |
114 |
117 /** |
115 /** |
118 * Checks an {@code int} upper bound value for validity. |
116 * Checks an {@code int} upper bound value for validity. |
119 * |
117 * |
120 * @param bound the upper bound (exclusive) |
118 * @param bound the upper bound (exclusive) |
|
119 * |
121 * @throws IllegalArgumentException if {@code bound} is not positive |
120 * @throws IllegalArgumentException if {@code bound} is not positive |
122 */ |
121 */ |
123 public static void checkBound(int bound) { |
122 public static void checkBound(int bound) { |
124 if (bound <= 0) |
123 if (bound <= 0) { |
125 throw new IllegalArgumentException(BadBound); |
124 throw new IllegalArgumentException(BAD_BOUND); |
|
125 } |
126 } |
126 } |
127 |
127 |
128 /** |
128 /** |
129 * Checks a {@code long} upper bound value for validity. |
129 * Checks a {@code long} upper bound value for validity. |
130 * |
130 * |
131 * @param bound the upper bound (exclusive) |
131 * @param bound the upper bound (exclusive) |
|
132 * |
132 * @throws IllegalArgumentException if {@code bound} is not positive |
133 * @throws IllegalArgumentException if {@code bound} is not positive |
133 */ |
134 */ |
134 public static void checkBound(long bound) { |
135 public static void checkBound(long bound) { |
135 if (bound <= 0) |
136 if (bound <= 0) { |
136 throw new IllegalArgumentException(BadBound); |
137 throw new IllegalArgumentException(BAD_BOUND); |
|
138 } |
137 } |
139 } |
138 |
140 |
139 /** |
141 /** |
140 * Checks a {@code float} range for validity. |
142 * Checks a {@code float} range for validity. |
141 * |
143 * |
142 * @param origin the least value (inclusive) in the range |
144 * @param origin the least value (inclusive) in the range |
143 * @param bound the upper bound (exclusive) of the range |
145 * @param bound the upper bound (exclusive) of the range |
144 * @throws IllegalArgumentException unless {@code origin} is finite, |
146 * |
145 * {@code bound} is finite, and {@code bound - origin} is finite |
147 * @throws IllegalArgumentException unless {@code origin} is finite, {@code bound} is finite, |
|
148 * and {@code bound - origin} is finite |
146 */ |
149 */ |
147 public static void checkRange(float origin, float bound) { |
150 public static void checkRange(float origin, float bound) { |
148 if (!(origin < bound && (bound - origin) < Float.POSITIVE_INFINITY)) |
151 if (!(origin < bound && (bound - origin) < Float.POSITIVE_INFINITY)) { |
149 throw new IllegalArgumentException(BadRange); |
152 throw new IllegalArgumentException(BAD_RANGE); |
|
153 } |
150 } |
154 } |
151 |
155 |
152 /** |
156 /** |
153 * Checks a {@code double} range for validity. |
157 * Checks a {@code double} range for validity. |
154 * |
158 * |
155 * @param origin the least value (inclusive) in the range |
159 * @param origin the least value (inclusive) in the range |
156 * @param bound the upper bound (exclusive) of the range |
160 * @param bound the upper bound (exclusive) of the range |
157 * @throws IllegalArgumentException unless {@code origin} is finite, |
161 * |
158 * {@code bound} is finite, and {@code bound - origin} is finite |
162 * @throws IllegalArgumentException unless {@code origin} is finite, {@code bound} is finite, |
|
163 * and {@code bound - origin} is finite |
159 */ |
164 */ |
160 public static void checkRange(double origin, double bound) { |
165 public static void checkRange(double origin, double bound) { |
161 if (!(origin < bound && (bound - origin) < Double.POSITIVE_INFINITY)) |
166 if (!(origin < bound && (bound - origin) < Double.POSITIVE_INFINITY)) { |
162 throw new IllegalArgumentException(BadRange); |
167 throw new IllegalArgumentException(BAD_RANGE); |
|
168 } |
163 } |
169 } |
164 |
170 |
165 /** |
171 /** |
166 * Checks an {@code int} range for validity. |
172 * Checks an {@code int} range for validity. |
167 * |
173 * |
168 * @param origin the least value that can be returned |
174 * @param origin the least value that can be returned |
169 * @param bound the upper bound (exclusive) for the returned value |
175 * @param bound the upper bound (exclusive) for the returned value |
170 * @throws IllegalArgumentException if {@code origin} is greater than |
176 * |
171 * or equal to {@code bound} |
177 * @throws IllegalArgumentException if {@code origin} is greater than or equal to {@code bound} |
172 */ |
178 */ |
173 public static void checkRange(int origin, int bound) { |
179 public static void checkRange(int origin, int bound) { |
174 if (origin >= bound) |
180 if (origin >= bound) { |
175 throw new IllegalArgumentException(BadRange); |
181 throw new IllegalArgumentException(BAD_RANGE); |
|
182 } |
176 } |
183 } |
177 |
184 |
178 /** |
185 /** |
179 * Checks a {@code long} range for validity. |
186 * Checks a {@code long} range for validity. |
180 * |
187 * |
181 * @param origin the least value that can be returned |
188 * @param origin the least value that can be returned |
182 * @param bound the upper bound (exclusive) for the returned value |
189 * @param bound the upper bound (exclusive) for the returned value |
183 * @throws IllegalArgumentException if {@code origin} is greater than |
190 * |
184 * or equal to {@code bound} |
191 * @throws IllegalArgumentException if {@code origin} is greater than or equal to {@code bound} |
185 */ |
192 */ |
186 public static void checkRange(long origin, long bound) { |
193 public static void checkRange(long origin, long bound) { |
187 if (origin >= bound) |
194 if (origin >= bound) { |
188 throw new IllegalArgumentException(BadRange); |
195 throw new IllegalArgumentException(BAD_RANGE); |
|
196 } |
189 } |
197 } |
190 |
198 |
191 /** |
199 /** |
192 * Given an array of seed bytes of any length, construct an array |
200 * Given an array of seed bytes of any length, construct an array |
193 * of {@code long} seed values of length {@code n}, such that the |
201 * of {@code long} seed values of length {@code n}, such that the |
237 * @param seed an array of {@code byte} values |
246 * @param seed an array of {@code byte} values |
238 * @param n the length of the result array (should be nonnegative) |
247 * @param n the length of the result array (should be nonnegative) |
239 * @param z the number of trailing result elements that are required |
248 * @param z the number of trailing result elements that are required |
240 * to be not all zero (should be nonnegative but not larger |
249 * to be not all zero (should be nonnegative but not larger |
241 * than {@code n}) |
250 * than {@code n}) |
|
251 * |
242 * @return an array of length {@code n} containing {@code int} seed values |
252 * @return an array of length {@code n} containing {@code int} seed values |
243 */ |
253 */ |
244 public static int[] convertSeedBytesToInts(byte[] seed, int n, int z) { |
254 public static int[] convertSeedBytesToInts(byte[] seed, int n, int z) { |
245 final int[] result = new int[n]; |
255 final int[] result = new int[n]; |
246 final int m = Math.min(seed.length, n << 2); |
256 final int m = Math.min(seed.length, n << 2); |
247 // Distribute seed bytes into the words to be formed. |
257 // Distribute seed bytes into the words to be formed. |
248 for (int j = 0; j < m; j++) { |
258 for (int j = 0; j < m; j++) { |
249 result[j>>2] = (result[j>>2] << 8) | seed[j]; |
259 result[j>>2] = (result[j>>2] << 8) | seed[j]; |
250 } |
260 } |
251 // If there aren't enough seed bytes for all the words we need, |
261 // If there aren't enough seed bytes for all the words we need, |
252 // use a SplitMix-style PRNG to fill in the rest. |
262 // use a SplitMix-style PRNG to fill in the rest. |
253 int v = result[0]; |
263 int v = result[0]; |
254 for (int j = (m + 3) >> 2; j < n; j++) { |
264 for (int j = (m + 3) >> 2; j < n; j++) { |
255 result[j] = mixMurmur32(v += SILVER_RATIO_32); |
265 result[j] = mixMurmur32(v += SILVER_RATIO_32); |
256 } |
266 } |
257 // Finally, we need to make sure the last z words are not all zero. |
267 // Finally, we need to make sure the last z words are not all zero. |
258 search: { |
268 search: { |
259 for (int j = n - z; j < n; j++) { |
269 for (int j = n - z; j < n; j++) { |
260 if (result[j] != 0) break search; |
270 if (result[j] != 0) break search; |
261 } |
271 } |
262 // If they are, fill in using a SplitMix-style PRNG. |
272 // If they are, fill in using a SplitMix-style PRNG. |
263 // Using "& ~1" in the next line defends against the case z==1 |
273 // Using "& ~1" in the next line defends against the case z==1 |
264 // by guaranteeing that the first generated value will be nonzero. |
274 // by guaranteeing that the first generated value will be nonzero. |
265 int w = result[0] & ~1; |
275 int w = result[0] & ~1; |
266 for (int j = n - z; j < n; j++) { |
276 for (int j = n - z; j < n; j++) { |
267 result[j] = mixMurmur32(w += SILVER_RATIO_32); |
277 result[j] = mixMurmur32(w += SILVER_RATIO_32); |
268 } |
278 } |
269 } |
279 } |
270 return result; |
280 return result; |
271 } |
281 } |
272 |
282 |
273 /* |
283 /* |
274 * Bounded versions of nextX methods used by streams, as well as |
284 * Bounded versions of nextX methods used by streams, as well as |
275 * the public nextX(origin, bound) methods. These exist mainly to |
285 * the public nextX(origin, bound) methods. These exist mainly to |
276 * avoid the need for multiple versions of stream spliterators |
286 * avoid the need for multiple versions of stream spliterators |
277 * across the different exported forms of streams. |
287 * across the different exported forms of streams. |
278 */ |
288 */ |
279 |
289 |
280 /** |
290 /** |
281 * This is the form of {@code nextLong} used by a {@code LongStream} |
291 * This is the form of {@code nextLong} used by a {@link LongStream} |
282 * {@code Spliterator} and by the public method |
292 * {@link Spliterator} and by the public method |
283 * {@code nextLong(origin, bound)}. If {@code origin} is greater |
293 * {@code nextLong(origin, bound)}. If {@code origin} is greater |
284 * than {@code bound}, then this method simply calls the unbounded |
294 * than {@code bound}, then this method simply calls the unbounded |
285 * version of {@code nextLong()}, choosing pseudorandomly from |
295 * version of {@code nextLong()}, choosing pseudorandomly from |
286 * among all 2<sup>64</sup> possible {@code long} values}, and |
296 * among all 2<sup>64</sup> possible {@code long} values}, and |
287 * otherwise uses one or more calls to {@code nextLong()} to |
297 * otherwise uses one or more calls to {@code nextLong()} to |
409 * </ol> |
420 * </ol> |
410 * |
421 * |
411 * @param rng a random number generator to be used as a |
422 * @param rng a random number generator to be used as a |
412 * source of pseudorandom {@code long} values |
423 * source of pseudorandom {@code long} values |
413 * @param bound the upper bound (exclusive); must be greater than zero |
424 * @param bound the upper bound (exclusive); must be greater than zero |
|
425 * |
414 * @return a pseudorandomly chosen {@code long} value |
426 * @return a pseudorandomly chosen {@code long} value |
415 */ |
427 */ |
416 public static long boundedNextLong(Rng rng, long bound) { |
428 public static long boundedNextLong(RandomNumberGenerator rng, long bound) { |
417 // Specialize boundedNextLong for origin == 0, bound > 0 |
429 // Specialize boundedNextLong for origin == 0, bound > 0 |
418 final long m = bound - 1; |
430 final long m = bound - 1; |
419 long r = rng.nextLong(); |
431 long r = rng.nextLong(); |
420 if ((bound & m) == 0L) { |
432 if ((bound & m) == 0L) { |
421 // The bound is a power of 2. |
433 // The bound is a power of 2. |
422 r &= m; |
434 r &= m; |
423 } else { |
435 } else { |
424 // Must reject over-represented candidates |
436 // Must reject over-represented candidates |
425 /* This loop takes an unlovable form (but it works): |
437 /* This loop takes an unlovable form (but it works): |
426 because the first candidate is already available, |
438 because the first candidate is already available, |
427 we need a break-in-the-middle construction, |
439 we need a break-in-the-middle construction, |
428 which is concisely but cryptically performed |
440 which is concisely but cryptically performed |
429 within the while-condition of a body-less for loop. */ |
441 within the while-condition of a body-less for loop. */ |
430 for (long u = r >>> 1; |
442 for (long u = r >>> 1; |
431 u + m - (r = u % bound) < 0L; |
443 u + m - (r = u % bound) < 0L; |
432 u = rng.nextLong() >>> 1) |
444 u = rng.nextLong() >>> 1) |
433 ; |
445 ; |
434 } |
446 } |
435 return r; |
447 return r; |
436 } |
448 } |
437 |
449 |
438 /** |
450 /** |
439 * This is the form of {@code nextInt} used by an {@code IntStream} |
451 * This is the form of {@code nextInt} used by an {@link IntStream} |
440 * {@code Spliterator} and by the public method |
452 * {@link Spliterator} and by the public method |
441 * {@code nextInt(origin, bound)}. If {@code origin} is greater |
453 * {@code nextInt(origin, bound)}. If {@code origin} is greater |
442 * than {@code bound}, then this method simply calls the unbounded |
454 * than {@code bound}, then this method simply calls the unbounded |
443 * version of {@code nextInt()}, choosing pseudorandomly from |
455 * version of {@code nextInt()}, choosing pseudorandomly from |
444 * among all 2<sup>64</sup> possible {@code int} values}, and |
456 * among all 2<sup>64</sup> possible {@code int} values}, and |
445 * otherwise uses one or more calls to {@code nextInt()} to |
457 * otherwise uses one or more calls to {@code nextInt()} to |
446 * choose a value pseudorandomly from the possible values |
458 * choose a value pseudorandomly from the possible values |
447 * between {@code origin} (inclusive) and {@code bound} (exclusive). |
459 * between {@code origin} (inclusive) and {@code bound} (exclusive). |
448 * |
460 * |
449 * @implNote The implementation of this method is identical to |
|
450 * the implementation of {@code nextLong(origin, bound)} |
|
451 * except that {@code int} values and the {@code nextInt()} |
|
452 * method are used rather than {@code long} values and the |
|
453 * {@code nextLong()} method. |
|
454 * |
|
455 * @param rng a random number generator to be used as a |
461 * @param rng a random number generator to be used as a |
456 * source of pseudorandom {@code int} values |
462 * source of pseudorandom {@code int} values |
457 * @param origin the least value that can be produced, |
463 * @param origin the least value that can be produced, |
458 * unless greater than or equal to {@code bound} |
464 * unless greater than or equal to {@code bound} |
459 * @param bound the upper bound (exclusive), unless {@code origin} |
465 * @param bound the upper bound (exclusive), unless {@code origin} |
460 * is greater than or equal to {@code bound} |
466 * is greater than or equal to {@code bound} |
|
467 * |
461 * @return a pseudorandomly chosen {@code int} value, |
468 * @return a pseudorandomly chosen {@code int} value, |
462 * which will be between {@code origin} (inclusive) and |
469 * which will be between {@code origin} (inclusive) and |
463 * {@code bound} exclusive unless {@code origin} |
470 * {@code bound} exclusive unless {@code origin} |
464 * is greater than or equal to {@code bound} |
471 * is greater than or equal to {@code bound} |
465 */ |
472 * |
466 public static int boundedNextInt(Rng rng, int origin, int bound) { |
473 * @implNote The implementation of this method is identical to |
|
474 * the implementation of {@code nextLong(origin, bound)} |
|
475 * except that {@code int} values and the {@code nextInt()} |
|
476 * method are used rather than {@code long} values and the |
|
477 * {@code nextLong()} method. |
|
478 */ |
|
479 public static int boundedNextInt(RandomNumberGenerator rng, int origin, int bound) { |
467 int r = rng.nextInt(); |
480 int r = rng.nextInt(); |
468 if (origin < bound) { |
481 if (origin < bound) { |
469 // It's not case (1). |
482 // It's not case (1). |
470 final int n = bound - origin; |
483 final int n = bound - origin; |
471 final int m = n - 1; |
484 final int m = n - 1; |
472 if ((n & m) == 0) { |
485 if ((n & m) == 0) { |
473 // It is case (2): length of range is a power of 2. |
486 // It is case (2): length of range is a power of 2. |
474 r = (r & m) + origin; |
487 r = (r & m) + origin; |
475 } else if (n > 0) { |
488 } else if (n > 0) { |
476 // It is case (3): need to reject over-represented candidates. |
489 // It is case (3): need to reject over-represented candidates. |
477 for (int u = r >>> 1; |
490 for (int u = r >>> 1; |
478 u + m - (r = u % n) < 0; |
491 u + m - (r = u % n) < 0; |
479 u = rng.nextInt() >>> 1) |
492 u = rng.nextInt() >>> 1) |
480 ; |
493 ; |
481 r += origin; |
494 r += origin; |
482 } |
495 } |
483 else { |
496 else { |
484 // It is case (4): length of range not representable as long. |
497 // It is case (4): length of range not representable as long. |
485 while (r < origin || r >= bound) |
498 while (r < origin || r >= bound) { |
486 |
499 r = rng.nextInt(); |
487 |
500 } |
488 r = rng.nextInt(); |
|
489 } |
501 } |
490 } |
502 } |
491 return r; |
503 return r; |
492 } |
504 } |
493 |
505 |
870 * low-order output bits do not have good statistical quality. |
883 * low-order output bits do not have good statistical quality. |
871 */ |
884 */ |
872 |
885 |
873 // Implementation support for nextExponential() |
886 // Implementation support for nextExponential() |
874 |
887 |
875 static double computeNextExponential(Rng rng) { |
888 static double computeNextExponential(RandomNumberGenerator rng) { |
876 long U1 = rng.nextLong(); |
889 long U1 = rng.nextLong(); |
877 // Experimentation on a variety of machines indicates that it is overall much faster |
890 // Experimentation on a variety of machines indicates that it is overall much faster |
878 // to do the following & and < operations on longs rather than first cast U1 to int |
891 // to do the following & and < operations on longs rather than first cast U1 to int |
879 // (but then we need to cast to int before doing the array indexing operation). |
892 // (but then we need to cast to int before doing the array indexing operation). |
880 long i = U1 & DoubleZigguratTables.exponentialLayerMask; |
893 long i = U1 & DoubleZigguratTables.exponentialLayerMask; |
881 if (i < DoubleZigguratTables.exponentialNumberOfLayers) { |
894 if (i < DoubleZigguratTables.exponentialNumberOfLayers) { |
882 // This is the fast path (occurring more than 98% of the time). Make an early exit. |
895 // This is the fast path (occurring more than 98% of the time). Make an early exit. |
883 return DoubleZigguratTables.exponentialX[(int)i] * (U1 >>> 1); |
896 return DoubleZigguratTables.exponentialX[(int)i] * (U1 >>> 1); |
884 } |
897 } |
885 // We didn't use the upper part of U1 after all. We'll be able to use it later. |
898 // We didn't use the upper part of U1 after all. We'll be able to use it later. |
886 |
899 |
887 for (double extra = 0.0; ; ) { |
900 for (double extra = 0.0; ; ) { |
888 // Use Walker's alias method to sample an (unsigned) integer j from a discrete |
901 // Use Walker's alias method to sample an (unsigned) integer j from a discrete |
889 // probability distribution that includes the tail and all the ziggurat overhangs; |
902 // probability distribution that includes the tail and all the ziggurat overhangs; |
890 // j will be less than DoubleZigguratTables.exponentialNumberOfLayers + 1. |
903 // j will be less than DoubleZigguratTables.exponentialNumberOfLayers + 1. |
891 long UA = rng.nextLong(); |
904 long UA = rng.nextLong(); |
892 int j = (int)UA & DoubleZigguratTables.exponentialAliasMask; |
905 int j = (int)UA & DoubleZigguratTables.exponentialAliasMask; |
893 if (UA >= DoubleZigguratTables.exponentialAliasThreshold[j]) { |
906 if (UA >= DoubleZigguratTables.exponentialAliasThreshold[j]) { |
894 j = DoubleZigguratTables.exponentialAliasMap[j] & DoubleZigguratTables.exponentialSignCorrectionMask; |
907 j = DoubleZigguratTables.exponentialAliasMap[j] & DoubleZigguratTables.exponentialSignCorrectionMask; |
895 } |
908 } |
896 if (j > 0) { // Sample overhang j |
909 if (j > 0) { // Sample overhang j |
897 // For the exponential distribution, every overhang is convex. |
910 // For the exponential distribution, every overhang is convex. |
898 final double[] X = DoubleZigguratTables.exponentialX; |
911 final double[] X = DoubleZigguratTables.exponentialX; |
899 final double[] Y = DoubleZigguratTables.exponentialY; |
912 final double[] Y = DoubleZigguratTables.exponentialY; |
900 for (;; U1 = (rng.nextLong() >>> 1)) { |
913 for (;; U1 = (rng.nextLong() >>> 1)) { |
901 long U2 = (rng.nextLong() >>> 1); |
914 long U2 = (rng.nextLong() >>> 1); |
902 // Compute the actual x-coordinate of the randomly chosen point. |
915 // Compute the actual x-coordinate of the randomly chosen point. |
903 double x = (X[j] * 0x1.0p63) + ((X[j-1] - X[j]) * (double)U1); |
916 double x = (X[j] * 0x1.0p63) + ((X[j-1] - X[j]) * (double)U1); |
904 // Does the point lie below the curve? |
917 // Does the point lie below the curve? |
905 long Udiff = U2 - U1; |
918 long Udiff = U2 - U1; |
906 if (Udiff < 0) { |
919 if (Udiff < 0) { |
907 // We picked a point in the upper-right triangle. None of those can be accepted. |
920 // We picked a point in the upper-right triangle. None of those can be accepted. |
908 // So remap the point into the lower-left triangle and try that. |
921 // So remap the point into the lower-left triangle and try that. |
909 // In effect, we swap U1 and U2, and invert the sign of Udiff. |
922 // In effect, we swap U1 and U2, and invert the sign of Udiff. |
910 Udiff = -Udiff; |
923 Udiff = -Udiff; |
911 U2 = U1; |
924 U2 = U1; |
912 U1 -= Udiff; |
925 U1 -= Udiff; |
913 } |
926 } |
914 if (Udiff >= DoubleZigguratTables.exponentialConvexMargin) { |
927 if (Udiff >= DoubleZigguratTables.exponentialConvexMargin) { |
915 return x + extra; // The chosen point is way below the curve; accept it. |
928 return x + extra; // The chosen point is way below the curve; accept it. |
916 } |
929 } |
917 // Compute the actual y-coordinate of the randomly chosen point. |
930 // Compute the actual y-coordinate of the randomly chosen point. |
918 double y = (Y[j] * 0x1.0p63) + ((Y[j] - Y[j-1]) * (double)U2); |
931 double y = (Y[j] * 0x1.0p63) + ((Y[j] - Y[j-1]) * (double)U2); |
919 // Now see how that y-coordinate compares to the curve |
932 // Now see how that y-coordinate compares to the curve |
920 if (y <= Math.exp(-x)) { |
933 if (y <= Math.exp(-x)) { |
921 return x + extra; // The chosen point is below the curve; accept it. |
934 return x + extra; // The chosen point is below the curve; accept it. |
922 } |
935 } |
923 // Otherwise, we reject this sample and have to try again. |
936 // Otherwise, we reject this sample and have to try again. |
924 } |
937 } |
925 } |
938 } |
926 // We are now committed to sampling from the tail. We could do a recursive call |
939 // We are now committed to sampling from the tail. We could do a recursive call |
927 // and then add X[0] but we save some time and stack space by using an iterative loop. |
940 // and then add X[0] but we save some time and stack space by using an iterative loop. |
928 extra += DoubleZigguratTables.exponentialX0; |
941 extra += DoubleZigguratTables.exponentialX0; |
929 // This is like the first five lines of this method, but if it returns, it first adds "extra". |
942 // This is like the first five lines of this method, but if it returns, it first adds "extra". |
930 U1 = rng.nextLong(); |
943 U1 = rng.nextLong(); |
931 i = U1 & DoubleZigguratTables.exponentialLayerMask; |
944 i = U1 & DoubleZigguratTables.exponentialLayerMask; |
932 if (i < DoubleZigguratTables.exponentialNumberOfLayers) { |
945 if (i < DoubleZigguratTables.exponentialNumberOfLayers) { |
933 return DoubleZigguratTables.exponentialX[(int)i] * (U1 >>> 1) + extra; |
946 return DoubleZigguratTables.exponentialX[(int)i] * (U1 >>> 1) + extra; |
934 } |
947 } |
935 } |
948 } |
936 } |
949 } |
937 |
950 |
938 // Implementation support for nextGaussian() |
951 // Implementation support for nextGaussian() |
939 |
952 |
940 static double computeNextGaussian(Rng rng) { |
953 static double computeNextGaussian(RandomNumberGenerator rng) { |
941 long U1 = rng.nextLong(); |
954 long U1 = rng.nextLong(); |
942 // Experimentation on a variety of machines indicates that it is overall much faster |
955 // Experimentation on a variety of machines indicates that it is overall much faster |
943 // to do the following & and < operations on longs rather than first cast U1 to int |
956 // to do the following & and < operations on longs rather than first cast U1 to int |
944 // (but then we need to cast to int before doing the array indexing operation). |
957 // (but then we need to cast to int before doing the array indexing operation). |
945 long i = U1 & DoubleZigguratTables.normalLayerMask; |
958 long i = U1 & DoubleZigguratTables.normalLayerMask; |
946 |
959 |
947 if (i < DoubleZigguratTables.normalNumberOfLayers) { |
960 if (i < DoubleZigguratTables.normalNumberOfLayers) { |
948 // This is the fast path (occurring more than 98% of the time). Make an early exit. |
961 // This is the fast path (occurring more than 98% of the time). Make an early exit. |
949 return DoubleZigguratTables.normalX[(int)i] * U1; // Note that the sign bit of U1 is used here. |
962 return DoubleZigguratTables.normalX[(int)i] * U1; // Note that the sign bit of U1 is used here. |
950 } |
963 } |
951 // We didn't use the upper part of U1 after all. |
964 // We didn't use the upper part of U1 after all. |
952 // Pull U1 apart into a sign bit and a 63-bit value for later use. |
965 // Pull U1 apart into a sign bit and a 63-bit value for later use. |
953 double signBit = (U1 >= 0) ? 1.0 : -1.0; |
966 double signBit = (U1 >= 0) ? 1.0 : -1.0; |
954 U1 = (U1 << 1) >>> 1; |
967 U1 = (U1 << 1) >>> 1; |
955 |
968 |
956 // Use Walker's alias method to sample an (unsigned) integer j from a discrete |
969 // Use Walker's alias method to sample an (unsigned) integer j from a discrete |
957 // probability distribution that includes the tail and all the ziggurat overhangs; |
970 // probability distribution that includes the tail and all the ziggurat overhangs; |
958 // j will be less than DoubleZigguratTables.normalNumberOfLayers + 1. |
971 // j will be less than DoubleZigguratTables.normalNumberOfLayers + 1. |
959 long UA = rng.nextLong(); |
972 long UA = rng.nextLong(); |
960 int j = (int)UA & DoubleZigguratTables.normalAliasMask; |
973 int j = (int)UA & DoubleZigguratTables.normalAliasMask; |
961 if (UA >= DoubleZigguratTables.normalAliasThreshold[j]) { |
974 if (UA >= DoubleZigguratTables.normalAliasThreshold[j]) { |
962 j = DoubleZigguratTables.normalAliasMap[j] & DoubleZigguratTables.normalSignCorrectionMask; |
975 j = DoubleZigguratTables.normalAliasMap[j] & DoubleZigguratTables.normalSignCorrectionMask; |
963 } |
976 } |
964 |
977 |
965 double x; |
978 double x; |
966 // Now the goal is to choose the result, which will be multiplied by signBit just before return. |
979 // Now the goal is to choose the result, which will be multiplied by signBit just before return. |
967 |
980 |
968 // There are four kinds of overhangs: |
981 // There are four kinds of overhangs: |
969 // |
982 // |
970 // j == 0 : Sample from tail |
983 // j == 0 : Sample from tail |
971 // 0 < j < normalInflectionIndex : Overhang is convex; can reject upper-right triangle |
984 // 0 < j < normalInflectionIndex : Overhang is convex; can reject upper-right triangle |
972 // j == normalInflectionIndex : Overhang includes the inflection point |
985 // j == normalInflectionIndex : Overhang includes the inflection point |
973 // j > normalInflectionIndex : Overhang is concave; can accept point in lower-left triangle |
986 // j > normalInflectionIndex : Overhang is concave; can accept point in lower-left triangle |
974 // |
987 // |
975 // Choose one of four loops to compute x, each specialized for a specific kind of overhang. |
988 // Choose one of four loops to compute x, each specialized for a specific kind of overhang. |
976 // Conditional statements are arranged such that the more likely outcomes are first. |
989 // Conditional statements are arranged such that the more likely outcomes are first. |
977 |
990 |
978 // In the three cases other than the tail case: |
991 // In the three cases other than the tail case: |
979 // U1 represents a fraction (scaled by 2**63) of the width of rectangle measured from the left. |
992 // U1 represents a fraction (scaled by 2**63) of the width of rectangle measured from the left. |
980 // U2 represents a fraction (scaled by 2**63) of the height of rectangle measured from the top. |
993 // U2 represents a fraction (scaled by 2**63) of the height of rectangle measured from the top. |
981 // Together they indicate a randomly chosen point within the rectangle. |
994 // Together they indicate a randomly chosen point within the rectangle. |
982 |
995 |
983 final double[] X = DoubleZigguratTables.normalX; |
996 final double[] X = DoubleZigguratTables.normalX; |
984 final double[] Y = DoubleZigguratTables.normalY; |
997 final double[] Y = DoubleZigguratTables.normalY; |
985 if (j > DoubleZigguratTables.normalInflectionIndex) { // Concave overhang |
998 if (j > DoubleZigguratTables.normalInflectionIndex) { // Concave overhang |
986 for (;; U1 = (rng.nextLong() >>> 1)) { |
999 for (;; U1 = (rng.nextLong() >>> 1)) { |
987 long U2 = (rng.nextLong() >>> 1); |
1000 long U2 = (rng.nextLong() >>> 1); |
988 // Compute the actual x-coordinate of the randomly chosen point. |
1001 // Compute the actual x-coordinate of the randomly chosen point. |
989 x = (X[j] * 0x1.0p63) + ((X[j-1] - X[j]) * (double)U1); |
1002 x = (X[j] * 0x1.0p63) + ((X[j-1] - X[j]) * (double)U1); |
990 // Does the point lie below the curve? |
1003 // Does the point lie below the curve? |
991 long Udiff = U2 - U1; |
1004 long Udiff = U2 - U1; |
992 if (Udiff >= 0) { |
1005 if (Udiff >= 0) { |
993 break; // The chosen point is in the lower-left triangle; accept it. |
1006 break; // The chosen point is in the lower-left triangle; accept it. |
994 } |
1007 } |
995 if (Udiff <= -DoubleZigguratTables.normalConcaveMargin) { |
1008 if (Udiff <= -DoubleZigguratTables.normalConcaveMargin) { |
996 continue; // The chosen point is way above the curve; reject it. |
1009 continue; // The chosen point is way above the curve; reject it. |
997 } |
1010 } |
998 // Compute the actual y-coordinate of the randomly chosen point. |
1011 // Compute the actual y-coordinate of the randomly chosen point. |
999 double y = (Y[j] * 0x1.0p63) + ((Y[j] - Y[j-1]) * (double)U2); |
1012 double y = (Y[j] * 0x1.0p63) + ((Y[j] - Y[j-1]) * (double)U2); |
1000 // Now see how that y-coordinate compares to the curve |
1013 // Now see how that y-coordinate compares to the curve |
1001 if (y <= Math.exp(-0.5*x*x)) { |
1014 if (y <= Math.exp(-0.5*x*x)) { |
1002 break; // The chosen point is below the curve; accept it. |
1015 break; // The chosen point is below the curve; accept it. |
1003 } |
1016 } |
1004 // Otherwise, we reject this sample and have to try again. |
1017 // Otherwise, we reject this sample and have to try again. |
1005 } |
1018 } |
1006 } else if (j == 0) { // Tail |
1019 } else if (j == 0) { // Tail |
1007 // Tail-sampling method of Marsaglia and Tsang. See any one of: |
1020 // Tail-sampling method of Marsaglia and Tsang. See any one of: |
1008 // Marsaglia and Tsang. 1984. A fast, easily implemented method for sampling from decreasing |
1021 // Marsaglia and Tsang. 1984. A fast, easily implemented method for sampling from decreasing |
1009 // or symmetric unimodal density functions. SIAM J. Sci. Stat. Comput. 5, 349-359. |
1022 // or symmetric unimodal density functions. SIAM J. Sci. Stat. Comput. 5, 349-359. |
1010 // Marsaglia and Tsang. 1998. The Monty Python method for generating random variables. |
1023 // Marsaglia and Tsang. 1998. The Monty Python method for generating random variables. |
1011 // ACM Trans. Math. Softw. 24, 3 (September 1998), 341-350. See page 342, step (4). |
1024 // ACM Trans. Math. Softw. 24, 3 (September 1998), 341-350. See page 342, step (4). |
1012 // http://doi.org/10.1145/292395.292453 |
1025 // http://doi.org/10.1145/292395.292453 |
1013 // Thomas, Luk, Leong, and Villasenor. 2007. Gaussian random number generators. |
1026 // Thomas, Luk, Leong, and Villasenor. 2007. Gaussian random number generators. |
1014 // ACM Comput. Surv. 39, 4, Article 11 (November 2007). See Algorithm 16. |
1027 // ACM Comput. Surv. 39, 4, Article 11 (November 2007). See Algorithm 16. |
1015 // http://doi.org/10.1145/1287620.1287622 |
1028 // http://doi.org/10.1145/1287620.1287622 |
1016 // Compute two separate random exponential samples and then compare them in certain way. |
1029 // Compute two separate random exponential samples and then compare them in certain way. |
1017 do { |
1030 do { |
1018 x = (1.0 / DoubleZigguratTables.normalX0) * computeNextExponential(rng); |
1031 x = (1.0 / DoubleZigguratTables.normalX0) * computeNextExponential(rng); |
1019 } while (computeNextExponential(rng) < 0.5*x*x); |
1032 } while (computeNextExponential(rng) < 0.5*x*x); |
1020 x += DoubleZigguratTables.normalX0; |
1033 x += DoubleZigguratTables.normalX0; |
1021 } else if (j < DoubleZigguratTables.normalInflectionIndex) { // Convex overhang |
1034 } else if (j < DoubleZigguratTables.normalInflectionIndex) { // Convex overhang |
1022 for (;; U1 = (rng.nextLong() >>> 1)) { |
1035 for (;; U1 = (rng.nextLong() >>> 1)) { |
1023 long U2 = (rng.nextLong() >>> 1); |
1036 long U2 = (rng.nextLong() >>> 1); |
1024 // Compute the actual x-coordinate of the randomly chosen point. |
1037 // Compute the actual x-coordinate of the randomly chosen point. |
1025 x = (X[j] * 0x1.0p63) + ((X[j-1] - X[j]) * (double)U1); |
1038 x = (X[j] * 0x1.0p63) + ((X[j-1] - X[j]) * (double)U1); |
1026 // Does the point lie below the curve? |
1039 // Does the point lie below the curve? |
1027 long Udiff = U2 - U1; |
1040 long Udiff = U2 - U1; |
1028 if (Udiff < 0) { |
1041 if (Udiff < 0) { |
1029 // We picked a point in the upper-right triangle. None of those can be accepted. |
1042 // We picked a point in the upper-right triangle. None of those can be accepted. |
1030 // So remap the point into the lower-left triangle and try that. |
1043 // So remap the point into the lower-left triangle and try that. |
1031 // In effect, we swap U1 and U2, and invert the sign of Udiff. |
1044 // In effect, we swap U1 and U2, and invert the sign of Udiff. |
1032 Udiff = -Udiff; |
1045 Udiff = -Udiff; |
1033 U2 = U1; |
1046 U2 = U1; |
1034 U1 -= Udiff; |
1047 U1 -= Udiff; |
1035 } |
1048 } |
1036 if (Udiff >= DoubleZigguratTables.normalConvexMargin) { |
1049 if (Udiff >= DoubleZigguratTables.normalConvexMargin) { |
1037 break; // The chosen point is way below the curve; accept it. |
1050 break; // The chosen point is way below the curve; accept it. |
1038 } |
1051 } |
1039 // Compute the actual y-coordinate of the randomly chosen point. |
1052 // Compute the actual y-coordinate of the randomly chosen point. |
1040 double y = (Y[j] * 0x1.0p63) + ((Y[j] - Y[j-1]) * (double)U2); |
1053 double y = (Y[j] * 0x1.0p63) + ((Y[j] - Y[j-1]) * (double)U2); |
1041 // Now see how that y-coordinate compares to the curve |
1054 // Now see how that y-coordinate compares to the curve |
1042 if (y <= Math.exp(-0.5*x*x)) break; // The chosen point is below the curve; accept it. |
1055 if (y <= Math.exp(-0.5*x*x)) break; // The chosen point is below the curve; accept it. |
1043 // Otherwise, we reject this sample and have to try again. |
1056 // Otherwise, we reject this sample and have to try again. |
1044 } |
1057 } |
1045 } else { |
1058 } else { |
1046 // The overhang includes the inflection point, so the curve is both convex and concave |
1059 // The overhang includes the inflection point, so the curve is both convex and concave |
1047 for (;; U1 = (rng.nextLong() >>> 1)) { |
1060 for (;; U1 = (rng.nextLong() >>> 1)) { |
1048 long U2 = (rng.nextLong() >>> 1); |
1061 long U2 = (rng.nextLong() >>> 1); |
1049 // Compute the actual x-coordinate of the randomly chosen point. |
1062 // Compute the actual x-coordinate of the randomly chosen point. |
1050 x = (X[j] * 0x1.0p63) + ((X[j-1] - X[j]) * (double)U1); |
1063 x = (X[j] * 0x1.0p63) + ((X[j-1] - X[j]) * (double)U1); |
1051 // Does the point lie below the curve? |
1064 // Does the point lie below the curve? |
1052 long Udiff = U2 - U1; |
1065 long Udiff = U2 - U1; |
1053 if (Udiff >= DoubleZigguratTables.normalConvexMargin) { |
1066 if (Udiff >= DoubleZigguratTables.normalConvexMargin) { |
1054 break; // The chosen point is way below the curve; accept it. |
1067 break; // The chosen point is way below the curve; accept it. |
1055 } |
1068 } |
1056 if (Udiff <= -DoubleZigguratTables.normalConcaveMargin) { |
1069 if (Udiff <= -DoubleZigguratTables.normalConcaveMargin) { |
1057 continue; // The chosen point is way above the curve; reject it. |
1070 continue; // The chosen point is way above the curve; reject it. |
1058 } |
1071 } |
1059 // Compute the actual y-coordinate of the randomly chosen point. |
1072 // Compute the actual y-coordinate of the randomly chosen point. |
1060 double y = (Y[j] * 0x1.0p63) + ((Y[j] - Y[j-1]) * (double)U2); |
1073 double y = (Y[j] * 0x1.0p63) + ((Y[j] - Y[j-1]) * (double)U2); |
1061 // Now see how that y-coordinate compares to the curve |
1074 // Now see how that y-coordinate compares to the curve |
1062 if (y <= Math.exp(-0.5*x*x)) { |
1075 if (y <= Math.exp(-0.5*x*x)) { |
1063 break; // The chosen point is below the curve; accept it. |
1076 break; // The chosen point is below the curve; accept it. |
1064 } |
1077 } |
1065 // Otherwise, we reject this sample and have to try again. |
1078 // Otherwise, we reject this sample and have to try again. |
1066 } |
1079 } |
1067 } |
1080 } |
1068 return signBit*x; |
1081 return signBit*x; |
1069 } |
1082 } |
1070 |
1083 |
1071 } |
1084 } |
1072 |
1085 |