src/java.base/share/classes/java/util/random/RNGSupport.java
branchJDK-8193209-branch
changeset 57437 f02ffcb61dce
parent 57436 b0c958c0e6c6
equal deleted inserted replaced
57436:b0c958c0e6c6 57437:f02ffcb61dce
    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
   196      * @param seed an array of {@code byte} values
   204      * @param seed an array of {@code byte} values
   197      * @param n the length of the result array (should be nonnegative)
   205      * @param n the length of the result array (should be nonnegative)
   198      * @param z the number of trailing result elements that are required
   206      * @param z the number of trailing result elements that are required
   199      *        to be not all zero (should be nonnegative but not larger
   207      *        to be not all zero (should be nonnegative but not larger
   200      *        than {@code n})
   208      *        than {@code n})
       
   209      *
   201      * @return an array of length {@code n} containing {@code long} seed values
   210      * @return an array of length {@code n} containing {@code long} seed values
   202      */
   211      */
   203     public static long[] convertSeedBytesToLongs(byte[] seed, int n, int z) {
   212     public static long[] convertSeedBytesToLongs(byte[] seed, int n, int z) {
   204 	final long[] result = new long[n];
   213         final long[] result = new long[n];
   205 	final int m = Math.min(seed.length, n << 3);
   214         final int m = Math.min(seed.length, n << 3);
   206 	// Distribute seed bytes into the words to be formed.
   215         // Distribute seed bytes into the words to be formed.
   207 	for (int j = 0; j < m; j++) {
   216         for (int j = 0; j < m; j++) {
   208 	    result[j>>3] = (result[j>>3] << 8) | seed[j];
   217             result[j>>3] = (result[j>>3] << 8) | seed[j];
   209 	}
   218         }
   210 	// If there aren't enough seed bytes for all the words we need,
   219         // If there aren't enough seed bytes for all the words we need,
   211 	// use a SplitMix-style PRNG to fill in the rest.
   220         // use a SplitMix-style PRNG to fill in the rest.
   212 	long v = result[0];
   221         long v = result[0];
   213 	for (int j = (m + 7) >> 3; j < n; j++) {
   222         for (int j = (m + 7) >> 3; j < n; j++) {
   214 	    result[j] = mixMurmur64(v += SILVER_RATIO_64);
   223             result[j] = mixMurmur64(v += SILVER_RATIO_64);
   215 	}
   224         }
   216 	// Finally, we need to make sure the last z words are not all zero.
   225         // Finally, we need to make sure the last z words are not all zero.
   217 	search: {
   226         search: {
   218 	    for (int j = n - z; j < n; j++) {
   227             for (int j = n - z; j < n; j++) {
   219 		if (result[j] != 0) break search;
   228                 if (result[j] != 0) break search;
   220 	    }
   229             }
   221 	    // If they are, fill in using a SplitMix-style PRNG.
   230             // If they are, fill in using a SplitMix-style PRNG.
   222 	    // Using "& ~1L" in the next line defends against the case z==1
   231             // Using "& ~1L" in the next line defends against the case z==1
   223 	    // by guaranteeing that the first generated value will be nonzero.
   232             // by guaranteeing that the first generated value will be nonzero.
   224 	    long w = result[0] & ~1L;
   233             long w = result[0] & ~1L;
   225 	    for (int j = n - z; j < n; j++) {
   234             for (int j = n - z; j < n; j++) {
   226 		result[j] = mixMurmur64(w += SILVER_RATIO_64);
   235                 result[j] = mixMurmur64(w += SILVER_RATIO_64);
   227 	    }
   236             }
   228 	}
   237         }
   229 	return result;
   238         return result;
   230     }
   239     }
   231 
   240 
   232     /**
   241     /**
   233      * Given an array of seed bytes of any length, construct an array
   242      * Given an array of seed bytes of any length, construct an array
   234      * of {@code int} seed values of length {@code n}, such that the
   243      * of {@code int} 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
   295      * possible long values is equally likely to be chosen).
   305      * possible long values is equally likely to be chosen).
   296      * Under some circumstances (when the specified range is not
   306      * Under some circumstances (when the specified range is not
   297      * a power of 2), {@code nextLong()} may be called additional times
   307      * a power of 2), {@code nextLong()} may be called additional times
   298      * to ensure that that the values in the specified range are
   308      * to ensure that that the values in the specified range are
   299      * equally likely to be chosen (provided the assumption holds).
   309      * equally likely to be chosen (provided the assumption holds).
   300      *
   310      * <p>
   301      * <p> The implementation considers four cases:
   311      * The implementation considers four cases:
   302      * <ol>
   312      * <ol>
   303      *
   313      *
   304      * <li> If the {@code} bound} is less than or equal to the {@code origin}
   314      * <li> If the {@code} bound} is less than or equal to the {@code origin}
   305      *      (indicated an unbounded form), the 64-bit {@code long} value
   315      *      (indicated an unbounded form), the 64-bit {@code long} value
   306      *      obtained from {@code nextLong()} is returned directly.
   316      *      obtained from {@code nextLong()} is returned directly.
   332      *        source of pseudorandom {@code long} values
   342      *        source of pseudorandom {@code long} values
   333      * @param origin the least value that can be produced,
   343      * @param origin the least value that can be produced,
   334      *        unless greater than or equal to {@code bound}
   344      *        unless greater than or equal to {@code bound}
   335      * @param bound the upper bound (exclusive), unless {@code origin}
   345      * @param bound the upper bound (exclusive), unless {@code origin}
   336      *        is greater than or equal to {@code bound}
   346      *        is greater than or equal to {@code bound}
       
   347      *
   337      * @return a pseudorandomly chosen {@code long} value,
   348      * @return a pseudorandomly chosen {@code long} value,
   338      *         which will be between {@code origin} (inclusive) and
   349      *         which will be between {@code origin} (inclusive) and
   339      *         {@code bound} exclusive unless {@code origin}
   350      *         {@code bound} exclusive unless {@code origin}
   340      *         is greater than or equal to {@code bound}
   351      *         is greater than or equal to {@code bound}
   341      */
   352      */
   342     public static long boundedNextLong(Rng rng, long origin, long bound) {
   353     public static long boundedNextLong(RandomNumberGenerator rng, long origin, long bound) {
   343         long r = rng.nextLong();
   354         long r = rng.nextLong();
   344         if (origin < bound) {
   355         if (origin < bound) {
   345 	    // It's not case (1).
   356             // It's not case (1).
   346             final long n = bound - origin;
   357             final long n = bound - origin;
   347 	    final long m = n - 1;
   358             final long m = n - 1;
   348             if ((n & m) == 0L) {
   359             if ((n & m) == 0L) {
   349 		// It is case (2): length of range is a power of 2.
   360                 // It is case (2): length of range is a power of 2.
   350                 r = (r & m) + origin;
   361                 r = (r & m) + origin;
   351 	    } else if (n > 0L) {
   362             } else if (n > 0L) {
   352 		// It is case (3): need to reject over-represented candidates.
   363                 // It is case (3): need to reject over-represented candidates.
   353 		/* This loop takes an unlovable form (but it works):
   364                 /* This loop takes an unlovable form (but it works):
   354 		   because the first candidate is already available,
   365                    because the first candidate is already available,
   355 		   we need a break-in-the-middle construction,
   366                    we need a break-in-the-middle construction,
   356 		   which is concisely but cryptically performed
   367                    which is concisely but cryptically performed
   357 		   within the while-condition of a body-less for loop. */
   368                    within the while-condition of a body-less for loop. */
   358                 for (long u = r >>> 1;            // ensure nonnegative
   369                 for (long u = r >>> 1;            // ensure nonnegative
   359                      u + m - (r = u % n) < 0L;    // rejection check
   370                      u + m - (r = u % n) < 0L;    // rejection check
   360                      u = rng.nextLong() >>> 1) // retry
   371                      u = rng.nextLong() >>> 1) // retry
   361                     ;
   372                     ;
   362                 r += origin;
   373                 r += origin;
   363             }
   374             }
   364             else {              
   375             else {
   365 		// It is case (4): length of range not representable as long.
   376                 // It is case (4): length of range not representable as long.
   366                 while (r < origin || r >= bound)
   377                 while (r < origin || r >= bound)
   367                     r = rng.nextLong();
   378                     r = rng.nextLong();
   368             }
   379             }
   369         }
   380         }
   370         return r;
   381         return r;
   386      * possible long values is equally likely to be chosen).
   397      * possible long values is equally likely to be chosen).
   387      * Under some circumstances (when the specified range is not
   398      * Under some circumstances (when the specified range is not
   388      * a power of 2), {@code nextLong()} may be called additional times
   399      * a power of 2), {@code nextLong()} may be called additional times
   389      * to ensure that that the values in the specified range are
   400      * to ensure that that the values in the specified range are
   390      * equally likely to be chosen (provided the assumption holds).
   401      * equally likely to be chosen (provided the assumption holds).
   391      *
   402      * <p>
   392      * <p> The implementation considers two cases:
   403      * The implementation considers two cases:
   393      * <ol>
   404      * <ol>
   394      *
   405      *
   395      * <li> If {@code bound} is an exact power of two 2<sup><i>m</i></sup>
   406      * <li> If {@code bound} is an exact power of two 2<sup><i>m</i></sup>
   396      *      for some integer <i>m</i>, then return the sum of {@code origin}
   407      *      for some integer <i>m</i>, then return the sum of {@code origin}
   397      *      and the <i>m</i> lowest-order bits of the value from
   408      *      and the <i>m</i> lowest-order bits of the value from
   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 
   498      * specialized for the case where the {@code origin} is zero
   510      * specialized for the case where the {@code origin} is zero
   499      * and the {@code bound} is greater than zero.  The value
   511      * and the {@code bound} is greater than zero.  The value
   500      * returned is chosen pseudorandomly from nonnegative integer
   512      * returned is chosen pseudorandomly from nonnegative integer
   501      * values less than {@code bound}.
   513      * values less than {@code bound}.
   502      *
   514      *
   503      * @implNote The implementation of this method is identical to
       
   504      *     the implementation of {@code nextLong(bound)}
       
   505      *     except that {@code int} values and the {@code nextInt()}
       
   506      *     method are used rather than {@code long} values and the
       
   507      *     {@code nextLong()} method.
       
   508      *
       
   509      * @param rng a random number generator to be used as a
   515      * @param rng a random number generator to be used as a
   510      *        source of pseudorandom {@code long} values
   516      *        source of pseudorandom {@code long} values
   511      * @param bound the upper bound (exclusive); must be greater than zero
   517      * @param bound the upper bound (exclusive); must be greater than zero
       
   518      *
   512      * @return a pseudorandomly chosen {@code long} value
   519      * @return a pseudorandomly chosen {@code long} value
   513      */
   520      *
   514     public static int boundedNextInt(Rng rng, int bound) {
   521      * @implNote The implementation of this method is identical to
       
   522      *           the implementation of {@code nextLong(bound)}
       
   523      *           except that {@code int} values and the {@code nextInt()}
       
   524      *           method are used rather than {@code long} values and the
       
   525      *           {@code nextLong()} method.
       
   526      */
       
   527     public static int boundedNextInt(RandomNumberGenerator rng, int bound) {
   515         // Specialize boundedNextInt for origin == 0, bound > 0
   528         // Specialize boundedNextInt for origin == 0, bound > 0
   516         final int m = bound - 1;
   529         final int m = bound - 1;
   517         int r = rng.nextInt();
   530         int r = rng.nextInt();
   518         if ((bound & m) == 0) {
   531         if ((bound & m) == 0) {
   519 	    // The bound is a power of 2.
   532             // The bound is a power of 2.
   520             r &= m;
   533             r &= m;
   521 	} else {
   534         } else {
   522 	    // Must reject over-represented candidates
   535             // Must reject over-represented candidates
   523             for (int u = r >>> 1;
   536             for (int u = r >>> 1;
   524                  u + m - (r = u % bound) < 0;
   537                  u + m - (r = u % bound) < 0;
   525                  u = rng.nextInt() >>> 1)
   538                  u = rng.nextInt() >>> 1)
   526                 ;
   539                 ;
   527         }
   540         }
   528         return r;
   541         return r;
   529     }
   542     }
   530     
   543 
   531     /**
   544     /**
   532      * This is the form of {@code nextDouble} used by a {@code DoubleStream}
   545      * This is the form of {@code nextDouble} used by a {@link DoubleStream}
   533      * {@code Spliterator} and by the public method
   546      * {@link Spliterator} and by the public method
   534      * {@code nextDouble(origin, bound)}.  If {@code origin} is greater
   547      * {@code nextDouble(origin, bound)}.  If {@code origin} is greater
   535      * than {@code bound}, then this method simply calls the unbounded
   548      * than {@code bound}, then this method simply calls the unbounded
   536      * version of {@code nextDouble()}, and otherwise scales and translates
   549      * version of {@code nextDouble()}, and otherwise scales and translates
   537      * the result of a call to {@code nextDouble()} so that it lies
   550      * the result of a call to {@code nextDouble()} so that it lies
   538      * between {@code origin} (inclusive) and {@code bound} (exclusive).
   551      * between {@code origin} (inclusive) and {@code bound} (exclusive).
   564      *         {@code bound} exclusive unless {@code origin}
   577      *         {@code bound} exclusive unless {@code origin}
   565      *         is greater than or equal to {@code bound},
   578      *         is greater than or equal to {@code bound},
   566      *         in which case it will be between 0.0 (inclusive)
   579      *         in which case it will be between 0.0 (inclusive)
   567      *         and 1.0 (exclusive)
   580      *         and 1.0 (exclusive)
   568      */
   581      */
   569     public static double boundedNextDouble(Rng rng, double origin, double bound) {
   582     public static double boundedNextDouble(RandomNumberGenerator rng, double origin, double bound) {
   570         double r = rng.nextDouble();
   583         double r = rng.nextDouble();
   571         if (origin < bound) {
   584         if (origin < bound) {
   572             r = r * (bound - origin) + origin;
   585             r = r * (bound - origin) + origin;
   573             if (r >= bound)  // may need to correct a rounding problem
   586             if (r >= bound)  // may need to correct a rounding problem
   574                 r = Double.longBitsToDouble(Double.doubleToLongBits(bound) - 1);
   587                 r = Double.longBitsToDouble(Double.doubleToLongBits(bound) - 1);
   575         }
   588         }
   576         return r;
   589         return r;
   577     }
   590     }
   578     
   591 
   579     /**
   592     /**
   580      * This is the form of {@code nextDouble} used by the public method
   593      * This is the form of {@code nextDouble} used by the public method
   581      * {@code nextDouble(bound)}.  This is essentially a version of
   594      * {@code nextDouble(bound)}.  This is essentially a version of
   582      * {@code boundedNextDouble(origin, bound)} that has been
   595      * {@code boundedNextDouble(origin, bound)} that has been
   583      * specialized for the case where the {@code origin} is zero
   596      * specialized for the case where the {@code origin} is zero
   594      * @param bound the upper bound (exclusive); must be finite and
   607      * @param bound the upper bound (exclusive); must be finite and
   595      *        greater than zero
   608      *        greater than zero
   596      * @return a pseudorandomly chosen {@code double} value
   609      * @return a pseudorandomly chosen {@code double} value
   597      *         between zero (inclusive) and {@code bound} (exclusive)
   610      *         between zero (inclusive) and {@code bound} (exclusive)
   598      */
   611      */
   599     public static double boundedNextDouble(Rng rng, double bound) {
   612     public static double boundedNextDouble(RandomNumberGenerator rng, double bound) {
   600         // Specialize boundedNextDouble for origin == 0, bound > 0
   613         // Specialize boundedNextDouble for origin == 0, bound > 0
   601         double r = rng.nextDouble();
   614         double r = rng.nextDouble();
   602 	r = r * bound;
   615         r = r * bound;
   603 	if (r >= bound)  // may need to correct a rounding problem
   616         if (r >= bound)  // may need to correct a rounding problem
   604 	    r = Double.longBitsToDouble(Double.doubleToLongBits(bound) - 1);
   617             r = Double.longBitsToDouble(Double.doubleToLongBits(bound) - 1);
   605         return r;
   618         return r;
   606     }
   619     }
   607 
   620 
   608     /**
   621     /**
   609      * This is the form of {@code nextFloat} used by a {@code FloatStream}
   622      * This is the form of {@code nextFloat} used by a {@code Stream<Float>}
   610      * {@code Spliterator} (if there were any) and by the public method
   623      * {@link Spliterator} (if there were any) and by the public method
   611      * {@code nextFloat(origin, bound)}.  If {@code origin} is greater
   624      * {@code nextFloat(origin, bound)}.  If {@code origin} is greater
   612      * than {@code bound}, then this method simply calls the unbounded
   625      * than {@code bound}, then this method simply calls the unbounded
   613      * version of {@code nextFloat()}, and otherwise scales and translates
   626      * version of {@code nextFloat()}, and otherwise scales and translates
   614      * the result of a call to {@code nextFloat()} so that it lies
   627      * the result of a call to {@code nextFloat()} so that it lies
   615      * between {@code origin} (inclusive) and {@code bound} (exclusive).
   628      * between {@code origin} (inclusive) and {@code bound} (exclusive).
   631      *         {@code bound} exclusive unless {@code origin}
   644      *         {@code bound} exclusive unless {@code origin}
   632      *         is greater than or equal to {@code bound},
   645      *         is greater than or equal to {@code bound},
   633      *         in which case it will be between 0.0 (inclusive)
   646      *         in which case it will be between 0.0 (inclusive)
   634      *         and 1.0 (exclusive)
   647      *         and 1.0 (exclusive)
   635      */
   648      */
   636     public static float boundedNextFloat(Rng rng, float origin, float bound) {
   649     public static float boundedNextFloat(RandomNumberGenerator rng, float origin, float bound) {
   637         float r = rng.nextFloat();
   650         float r = rng.nextFloat();
   638         if (origin < bound) {
   651         if (origin < bound) {
   639             r = r * (bound - origin) + origin;
   652             r = r * (bound - origin) + origin;
   640             if (r >= bound) // may need to correct a rounding problem
   653             if (r >= bound) // may need to correct a rounding problem
   641                 r = Float.intBitsToFloat(Float.floatToIntBits(bound) - 1);
   654                 r = Float.intBitsToFloat(Float.floatToIntBits(bound) - 1);
   661      * @param bound the upper bound (exclusive); must be finite and
   674      * @param bound the upper bound (exclusive); must be finite and
   662      *        greater than zero
   675      *        greater than zero
   663      * @return a pseudorandomly chosen {@code float} value
   676      * @return a pseudorandomly chosen {@code float} value
   664      *         between zero (inclusive) and {@code bound} (exclusive)
   677      *         between zero (inclusive) and {@code bound} (exclusive)
   665      */
   678      */
   666     public static float boundedNextFloat(Rng rng, float bound) {
   679     public static float boundedNextFloat(RandomNumberGenerator rng, float bound) {
   667         // Specialize boundedNextFloat for origin == 0, bound > 0
   680         // Specialize boundedNextFloat for origin == 0, bound > 0
   668         float r = rng.nextFloat();
   681         float r = rng.nextFloat();
   669 	r = r * bound;
   682         r = r * bound;
   670 	if (r >= bound) // may need to correct a rounding problem
   683         if (r >= bound) // may need to correct a rounding problem
   671 	    r = Float.intBitsToFloat(Float.floatToIntBits(bound) - 1);
   684             r = Float.intBitsToFloat(Float.floatToIntBits(bound) - 1);
   672         return r;
   685         return r;
   673     }
   686     }
   674 
   687 
   675     // The following decides which of two strategies initialSeed() will use.
   688     // The following decides which of two strategies initialSeed() will use.
   676     private static boolean secureRandomSeedRequested() {
   689     private static boolean secureRandomSeedRequested() {
   677 	String pp = java.security.AccessController.doPrivileged(
   690         String pp = java.security.AccessController.doPrivileged(
   678                 new sun.security.action.GetPropertyAction(
   691                 new sun.security.action.GetPropertyAction(
   679                         "java.util.secureRandomSeed"));
   692                         "java.util.secureRandomSeed"));
   680 	return (pp != null && pp.equalsIgnoreCase("true"));
   693         return (pp != null && pp.equalsIgnoreCase("true"));
   681     }
   694     }
   682     
   695 
   683     private static final boolean useSecureRandomSeed = secureRandomSeedRequested();
   696     private static final boolean useSecureRandomSeed = secureRandomSeedRequested();
   684 	
   697 
   685     /**
   698     /**
   686      * Returns a {@code long} value (chosen from some
   699      * Returns a {@code long} value (chosen from some
   687      * machine-dependent entropy source) that may be useful for
   700      * machine-dependent entropy source) that may be useful for
   688      * initializing a source of seed values for instances of {@code Rng}
   701      * initializing a source of seed values for instances of {@link RandomNumberGenerator}
   689      * created by zero-argument constructors.  (This method should
   702      * created by zero-argument constructors.  (This method should
   690      * <i>not</i> be called repeatedly, once per constructed
   703      * <i>not</i> be called repeatedly, once per constructed
   691      * object; at most it should be called once per class.)
   704      * object; at most it should be called once per class.)
   692      *
   705      *
   693      * @return a {@code long} value, randomly chosen using
   706      * @return a {@code long} value, randomly chosen using
   702             return s;
   715             return s;
   703         }
   716         }
   704         return (mixStafford13(System.currentTimeMillis()) ^
   717         return (mixStafford13(System.currentTimeMillis()) ^
   705                 mixStafford13(System.nanoTime()));
   718                 mixStafford13(System.nanoTime()));
   706     }
   719     }
   707     
   720 
   708     /**
   721     /**
   709      * The first 32 bits of the golden ratio (1+sqrt(5))/2, forced to be odd.
   722      * The first 32 bits of the golden ratio (1+sqrt(5))/2, forced to be odd.
   710      * Useful for producing good Weyl sequences or as an arbitrary nonzero odd value.
   723      * Useful for producing good Weyl sequences or as an arbitrary nonzero odd value.
   711      */
   724      */
   712     public static final int  GOLDEN_RATIO_32 = 0x9e3779b9;
   725     public static final int  GOLDEN_RATIO_32 = 0x9e3779b9;
   713     
   726 
   714     /**
   727     /**
   715      * The first 64 bits of the golden ratio (1+sqrt(5))/2, forced to be odd.
   728      * The first 64 bits of the golden ratio (1+sqrt(5))/2, forced to be odd.
   716      * Useful for producing good Weyl sequences or as an arbitrary nonzero odd value.
   729      * Useful for producing good Weyl sequences or as an arbitrary nonzero odd value.
   717      */
   730      */
   718     public static final long GOLDEN_RATIO_64 = 0x9e3779b97f4a7c15L;
   731     public static final long GOLDEN_RATIO_64 = 0x9e3779b97f4a7c15L;
   719     
   732 
   720     /**
   733     /**
   721      * The first 32 bits of the silver ratio 1+sqrt(2), forced to be odd.
   734      * The first 32 bits of the silver ratio 1+sqrt(2), forced to be odd.
   722      * Useful for producing good Weyl sequences or as an arbitrary nonzero odd value.
   735      * Useful for producing good Weyl sequences or as an arbitrary nonzero odd value.
   723      */
   736      */
   724     public static final int  SILVER_RATIO_32 = 0x6A09E667;
   737     public static final int  SILVER_RATIO_32 = 0x6A09E667;
   725     
   738 
   726     /**
   739     /**
   727      * The first 64 bits of the silver ratio 1+sqrt(2), forced to be odd.
   740      * The first 64 bits of the silver ratio 1+sqrt(2), forced to be odd.
   728      * Useful for producing good Weyl sequences or as an arbitrary nonzero odd value.
   741      * Useful for producing good Weyl sequences or as an arbitrary nonzero odd value.
   729      */
   742      */
   730     public static final long SILVER_RATIO_64 = 0x6A09E667F3BCC909L;
   743     public static final long SILVER_RATIO_64 = 0x6A09E667F3BCC909L;
   731     
   744 
   732     /**
   745     /**
   733      * Computes the 64-bit mixing function for MurmurHash3.
   746      * Computes the 64-bit mixing function for MurmurHash3.
   734      * This is a 64-bit hashing function with excellent avalanche statistics.
   747      * This is a 64-bit hashing function with excellent avalanche statistics.
   735      * https://github.com/aappleby/smhasher/wiki/MurmurHash3
   748      * https://github.com/aappleby/smhasher/wiki/MurmurHash3
   736      *
   749      *
   743     public static long mixMurmur64(long z) {
   756     public static long mixMurmur64(long z) {
   744         z = (z ^ (z >>> 33)) * 0xff51afd7ed558ccdL;
   757         z = (z ^ (z >>> 33)) * 0xff51afd7ed558ccdL;
   745         z = (z ^ (z >>> 33)) * 0xc4ceb9fe1a85ec53L;
   758         z = (z ^ (z >>> 33)) * 0xc4ceb9fe1a85ec53L;
   746         return z ^ (z >>> 33);
   759         return z ^ (z >>> 33);
   747     }
   760     }
   748     
   761 
   749     /**
   762     /**
   750      * Computes Stafford variant 13 of the 64-bit mixing function for MurmurHash3.
   763      * Computes Stafford variant 13 of the 64-bit mixing function for MurmurHash3.
   751      * This is a 64-bit hashing function with excellent avalanche statistics.
   764      * This is a 64-bit hashing function with excellent avalanche statistics.
   752      * http://zimbry.blogspot.com/2011/09/better-bit-mixing-improving-on.html
   765      * http://zimbry.blogspot.com/2011/09/better-bit-mixing-improving-on.html
   753      *
   766      *
   760     public static long mixStafford13(long z) {
   773     public static long mixStafford13(long z) {
   761         z = (z ^ (z >>> 30)) * 0xbf58476d1ce4e5b9L;
   774         z = (z ^ (z >>> 30)) * 0xbf58476d1ce4e5b9L;
   762         z = (z ^ (z >>> 27)) * 0x94d049bb133111ebL;
   775         z = (z ^ (z >>> 27)) * 0x94d049bb133111ebL;
   763         return z ^ (z >>> 31);
   776         return z ^ (z >>> 31);
   764     }
   777     }
   765     
   778 
   766     /**
   779     /**
   767      * Computes Doug Lea's 64-bit mixing function.
   780      * Computes Doug Lea's 64-bit mixing function.
   768      * This is a 64-bit hashing function with excellent avalanche statistics.
   781      * This is a 64-bit hashing function with excellent avalanche statistics.
   769      * It has the advantages of using the same multiplicative constant twice
   782      * It has the advantages of using the same multiplicative constant twice
   770      * and of using only 32-bit shifts.
   783      * and of using only 32-bit shifts.
   814         z = (z ^ (z >>> 16)) * 0xd36d884b;
   827         z = (z ^ (z >>> 16)) * 0xd36d884b;
   815         z = (z ^ (z >>> 16)) * 0xd36d884b;
   828         z = (z ^ (z >>> 16)) * 0xd36d884b;
   816         return z ^ (z >>> 16);
   829         return z ^ (z >>> 16);
   817     }
   830     }
   818 
   831 
   819     // Non-public (package only) support for spliterators needed by AbstractSplittableRng
   832     // Non-public (package only) support for spliterators needed by AbstractSplittableRNG
   820     // and AbstractArbitrarilyJumpableRng and AbstractSharedRng
   833     // and AbstractArbitrarilyJumpableRNG and AbstractSharedRNG
   821 
   834 
   822     /**
   835     /**
   823      * Base class for making Spliterator classes for streams of randomly chosen values.
   836      * Base class for making Spliterator classes for streams of randomly chosen values.
   824      */
   837      */
   825     static abstract class RandomSpliterator {
   838     static abstract class RandomSpliterator {
   826         long index;
   839         long index;
   827         final long fence;
   840         final long fence;
   828 
   841 
   829 	RandomSpliterator(long index, long fence) {
   842         RandomSpliterator(long index, long fence) {
   830 	    this.index = index; this.fence = fence;
   843             this.index = index; this.fence = fence;
   831         }
   844         }
   832 	
   845 
   833         public long estimateSize() {
   846         public long estimateSize() {
   834             return fence - index;
   847             return fence - index;
   835         }
   848         }
   836 
   849 
   837         public int characteristics() {
   850         public int characteristics() {
   838             return (Spliterator.SIZED | Spliterator.SUBSIZED |
   851             return (Spliterator.SIZED | Spliterator.SUBSIZED |
   839                     Spliterator.NONNULL | Spliterator.IMMUTABLE);
   852                     Spliterator.NONNULL | Spliterator.IMMUTABLE);
   840         }
   853         }
   841     }
   854     }
   842     
   855 
   843     
   856 
   844     /* 
   857     /*
   845      * Implementation support for nextExponential() and nextGaussian() methods of Rng.
   858      * Implementation support for nextExponential() and nextGaussian() methods of RandomNumberGenerator.
   846      *
   859      *
   847      * Each is implemented using McFarland's fast modified ziggurat algorithm (largely
   860      * Each is implemented using McFarland's fast modified ziggurat algorithm (largely
   848      * table-driven, with rare cases handled by computation and rejection sampling).
   861      * table-driven, with rare cases handled by computation and rejection sampling).
   849      * Walker's alias method for sampling a discrete distribution also plays a role.
   862      * Walker's alias method for sampling a discrete distribution also plays a role.
   850      *
   863      *
   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