8273056: java.util.random does not correctly sample exponential or Gaussian distributions

Co-authored-by: Guy Steele <gls@openjdk.org>
Reviewed-by: bpb, darcy
This commit is contained in:
Jim Laskey 2021-12-02 13:00:14 +00:00
parent b79554bb5c
commit 3d98ec1b7b

View file

@ -1186,10 +1186,10 @@ public class RandomSupport {
// For the exponential distribution, every overhang is convex. // For the exponential distribution, every overhang is convex.
final double[] X = DoubleZigguratTables.exponentialX; final double[] X = DoubleZigguratTables.exponentialX;
final double[] Y = DoubleZigguratTables.exponentialY; 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); 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? // Does the point lie below the curve?
long Udiff = U2 - U1; long Udiff = U2 - U1;
if (Udiff < 0) { if (Udiff < 0) {
@ -1200,11 +1200,13 @@ public class RandomSupport {
U2 = U1; U2 = U1;
U1 -= Udiff; 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) { if (Udiff >= DoubleZigguratTables.exponentialConvexMargin) {
return x + extra; // The chosen point is way below the curve; accept it. return x + extra; // The chosen point is way below the curve; accept it.
} }
// Compute the actual y-coordinate of the randomly chosen point. // 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 // Now see how that y-coordinate compares to the curve
if (y <= Math.exp(-x)) { if (y <= Math.exp(-x)) {
return x + extra; // The chosen point is below the curve; accept it. return x + extra; // The chosen point is below the curve; accept it.
@ -1323,7 +1325,7 @@ public class RandomSupport {
continue; // The chosen point is way above the curve; reject it. continue; // The chosen point is way above the curve; reject it.
} }
// Compute the actual y-coordinate of the randomly chosen point. // 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 // Now see how that y-coordinate compares to the curve
if (y <= Math.exp(-0.5*x*x)) { if (y <= Math.exp(-0.5*x*x)) {
break; // The chosen point is below the curve; accept it. break; // The chosen point is below the curve; accept it.
@ -1348,8 +1350,6 @@ public class RandomSupport {
} else if (j < DoubleZigguratTables.normalInflectionIndex) { // Convex overhang } else if (j < DoubleZigguratTables.normalInflectionIndex) { // Convex overhang
for (;; U1 = (rng.nextLong() >>> 1)) { for (;; U1 = (rng.nextLong() >>> 1)) {
long U2 = (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? // Does the point lie below the curve?
long Udiff = U2 - U1; long Udiff = U2 - U1;
if (Udiff < 0) { if (Udiff < 0) {
@ -1360,11 +1360,13 @@ public class RandomSupport {
U2 = U1; U2 = U1;
U1 -= Udiff; 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) { if (Udiff >= DoubleZigguratTables.normalConvexMargin) {
break; // The chosen point is way below the curve; accept it. break; // The chosen point is way below the curve; accept it.
} }
// Compute the actual y-coordinate of the randomly chosen point. // 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 // 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. 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. // Otherwise, we reject this sample and have to try again.
@ -1384,7 +1386,7 @@ public class RandomSupport {
continue; // The chosen point is way above the curve; reject it. continue; // The chosen point is way above the curve; reject it.
} }
// Compute the actual y-coordinate of the randomly chosen point. // 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 // Now see how that y-coordinate compares to the curve
if (y <= Math.exp(-0.5*x*x)) { if (y <= Math.exp(-0.5*x*x)) {
break; // The chosen point is below the curve; accept it. break; // The chosen point is below the curve; accept it.