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

Backport-of: 3d98ec1b7bc77237177ecfc069c0e9a7e75829bc

Author: Jim Laskey

src/java.base/share/classes/jdk/internal/util/random/RandomSupport.java

@@ -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) {
@@ -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.
@@ -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.
@@ -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) {
@@ -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.
@@ -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.