/ jdk Public

## Commit

Browse the repository at this point in the history
8273056: java.util.random does not correctly sample exponential or Ga…
```…ussian distributions

Co-authored-by: Guy Steele <gls@openjdk.org>
Reviewed-by: bpb, darcy```
JimLaskey and Guy Steele committed Dec 2, 2021
1 parent b79554b commit 3d98ec1
Showing 1 changed file with 11 additions and 9 deletions.
20 changes: 11 additions & 9 deletions
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
Expand Up @@ -1186,10 +1186,10 @@ public static double computeNextExponential(RandomGenerator rng) {
// For the exponential distribution, every overhang is convex.
final double[] X = DoubleZigguratTables.exponentialX;
final double[] Y = DoubleZigguratTables.exponentialY;
for (;; U1 = (rng.nextLong() >>> 1)) {
// At this point, the high-order bits of U1 have not been used yet,
// but we need the value in U1 to be positive.
for (U1 = (U1 >>> 1);; U1 = (rng.nextLong() >>> 1)) {
long U2 = (rng.nextLong() >>> 1);
// Compute the actual x-coordinate of the randomly chosen point.
double x = (X[j] * 0x1.0p63) + ((X[j-1] - X[j]) * (double)U1);
// Does the point lie below the curve?
long Udiff = U2 - U1;
if (Udiff < 0) {
Expand All @@ -1200,11 +1200,13 @@ public static double computeNextExponential(RandomGenerator rng) {
U2 = U1;
U1 -= Udiff;
}
// Compute the actual x-coordinate of the randomly chosen point.
double x = (X[j] * 0x1.0p63) + ((X[j-1] - X[j]) * (double)U1);
if (Udiff >= DoubleZigguratTables.exponentialConvexMargin) {
return x + extra; // The chosen point is way below the curve; accept it.
}
// Compute the actual y-coordinate of the randomly chosen point.
double y = (Y[j] * 0x1.0p63) + ((Y[j] - Y[j-1]) * (double)U2);
double y = (Y[j] * 0x1.0p63) + ((Y[j-1] - Y[j]) * (double)U2);
// Now see how that y-coordinate compares to the curve
if (y <= Math.exp(-x)) {
return x + extra; // The chosen point is below the curve; accept it.
Expand Down Expand Up @@ -1323,7 +1325,7 @@ public static double computeNextGaussian(RandomGenerator rng) {
continue; // The chosen point is way above the curve; reject it.
}
// Compute the actual y-coordinate of the randomly chosen point.
double y = (Y[j] * 0x1.0p63) + ((Y[j] - Y[j-1]) * (double)U2);
double y = (Y[j] * 0x1.0p63) + ((Y[j-1] - Y[j]) * (double)U2);
// Now see how that y-coordinate compares to the curve
if (y <= Math.exp(-0.5*x*x)) {
break; // The chosen point is below the curve; accept it.
Expand All @@ -1348,8 +1350,6 @@ public static double computeNextGaussian(RandomGenerator rng) {
} else if (j < DoubleZigguratTables.normalInflectionIndex) { // Convex overhang
for (;; U1 = (rng.nextLong() >>> 1)) {
long U2 = (rng.nextLong() >>> 1);
// Compute the actual x-coordinate of the randomly chosen point.
x = (X[j] * 0x1.0p63) + ((X[j-1] - X[j]) * (double)U1);
// Does the point lie below the curve?
long Udiff = U2 - U1;
if (Udiff < 0) {
Expand All @@ -1360,11 +1360,13 @@ public static double computeNextGaussian(RandomGenerator rng) {
U2 = U1;
U1 -= Udiff;
}
// Compute the actual x-coordinate of the randomly chosen point.
x = (X[j] * 0x1.0p63) + ((X[j-1] - X[j]) * (double)U1);
if (Udiff >= DoubleZigguratTables.normalConvexMargin) {
break; // The chosen point is way below the curve; accept it.
}
// Compute the actual y-coordinate of the randomly chosen point.
double y = (Y[j] * 0x1.0p63) + ((Y[j] - Y[j-1]) * (double)U2);
double y = (Y[j] * 0x1.0p63) + ((Y[j-1] - Y[j]) * (double)U2);
// Now see how that y-coordinate compares to the curve
if (y <= Math.exp(-0.5*x*x)) break; // The chosen point is below the curve; accept it.
// Otherwise, we reject this sample and have to try again.
Expand All @@ -1384,7 +1386,7 @@ public static double computeNextGaussian(RandomGenerator rng) {
continue; // The chosen point is way above the curve; reject it.
}
// Compute the actual y-coordinate of the randomly chosen point.
double y = (Y[j] * 0x1.0p63) + ((Y[j] - Y[j-1]) * (double)U2);
double y = (Y[j] * 0x1.0p63) + ((Y[j-1] - Y[j]) * (double)U2);
// Now see how that y-coordinate compares to the curve
if (y <= Math.exp(-0.5*x*x)) {
break; // The chosen point is below the curve; accept it.
Expand Down