8302987: Add uniform and spatially equidistributed bounded double streams to RandomGenerator

Reviewed-by: darcy
This commit is contained in:
Raffaello Giulietti 2023-07-19 16:48:54 +00:00
parent d1c788c52b
commit 14cf035681
2 changed files with 520 additions and 2 deletions

View File

@ -1,5 +1,5 @@
/*
* Copyright (c) 2021, 2022, Oracle and/or its affiliates. All rights reserved.
* Copyright (c) 2021, 2023, Oracle and/or its affiliates. All rights reserved.
* DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
*
* This code is free software; you can redistribute it and/or modify it
@ -30,13 +30,14 @@ import java.security.SecureRandom;
import java.util.Objects;
import java.util.concurrent.ThreadLocalRandom;
import jdk.internal.util.random.RandomSupport;
import jdk.internal.util.random.RandomSupport.*;
import java.util.stream.DoubleStream;
import java.util.stream.IntStream;
import java.util.stream.LongStream;
import java.util.stream.Stream;
import static java.lang.Math.*;
/**
* The {@link RandomGenerator} interface is designed to provide a common
* protocol for objects that generate random or (more typically) pseudorandom
@ -252,6 +253,174 @@ public interface RandomGenerator {
return doubles(randomNumberOrigin, randomNumberBound).limit(streamSize);
}
/**
* Returns an effectively unlimited stream of pseudorandomly chosen
* {@code double} values, where each value is between the specified
* {@code left} boundary and the specified {@code right} boundary.
* The {@code left} boundary is included as indicated by
* {@code isLeftIncluded}; similarly, the {@code right} boundary is included
* as indicated by {@code isRightIncluded}.
*
* <p>The stream potentially produces all multiples <i>k</i> &delta;
* (<i>k</i> integer) lying in the interval specified by the parameters,
* where &delta; > 0 is the smallest number for which all these multiples
* are exact {@code double}s.
* They are therefore all equidistant.
* The uniformity of the distribution of the {@code double}s produced by
* the stream depends on the quality of the underlying {@link #nextLong(long)}.
*
* @implSpec The default implementation first determines the &delta; above.
* It then computes both the smallest integer <i>k</i><sub><i>l</i></sub>
* such that <i>k</i><sub><i>l</i></sub> &delta; lies <em>inside</em>
* the given interval, and the smallest integer <i>n</i> > 0 such that
* (<i>k</i><sub><i>l</i></sub> + <i>n</i>) &delta; lies
* <em>outside</em> the interval.
* Finally, it returns a stream which generates the {@code double}s
* according to (<i>k</i><sub><i>l</i></sub> + {@code nextLong(}<i>n</i>{@code )})
* &delta;.
* The stream never produces {@code -0.0}, although it may produce
* {@code 0.0} if the specified interval contains 0.
*
* @param left the left boundary
* @param right the right boundary
* @param isLeftIncluded whether the {@code left} boundary is included
* @param isRightIncluded whether the {@code right} boundary is included
*
* @return a stream of pseudorandomly chosen {@code double} values, each
* between {@code left} and {@code right}, as specified above.
*
* @throws IllegalArgumentException if {@code left} is not finite,
* or {@code right} is not finite, or if the specified interval
* is empty.
*
* @since 21
*/
default DoubleStream equiDoubles(double left, double right,
boolean isLeftIncluded, boolean isRightIncluded) {
if (!(Double.NEGATIVE_INFINITY < left
&& right < Double.POSITIVE_INFINITY
&& (isLeftIncluded ? left : nextUp(left))
<= (isRightIncluded ? right : nextDown(right)))) {
throw new IllegalArgumentException(
"the boundaries must be finite and the interval must not be empty");
}
/*
* Inspired by
* Goualard, "Drawing random floating-point numbers from an
* interval", ACM TOMACS, 2022, 32 (3)
* (https://hal.science/hal-03282794v4)
* although implemented differently.
*
* It is assumed that left <= right.
* Whether the boundaries of the interval I = <left, right> are included
* is indicated by isLeftIncluded and isRightIncluded.
*
* delta > 0 is the smallest double such that every product k delta
* (k integer) that lies in I is an exact double as well.
* It turns out that delta is always a power of 2.
*
* kl is the smallest k such that k delta is inside I.
* kr > kl is the smallest k such that k delta is outside I.
* n is kr - kl
*/
double delta; // captured
long kl; // captured
long kr;
long n; // captured
if (left <= -right) {
/*
* Here,
* left <= 0, left <= right <= -left
* P = Double.PRECISION
*
* delta is the distance from left to the next double in the
* direction of positive infinity.
* Most of the time, this is equivalent to the ulp of left, but not
* always.
* For example, for left == -1.0, Math.ulp(left) == 2.220446049250313E-16,
* whereas delta as computed here is 1.1102230246251565E-16.
*
* Every product k delta lying in [left, -left] is an exact double.
* Thus, every product k delta lying in I is an exact double, too.
* Any other positive eps < delta does not meet this property:
* some product k eps lying in I is not an exact double.
* On the other hand, any other eps > delta would generate more
* sparse products k eps, that is, fewer doubles in I.
* delta is therefore the best value to ensure the largest number
* of equidistant doubles in the interval I.
*
* left / delta is an exact double and an exact integer with
* -2^P <= left / delta <= 0
* Thus, kl is computed exactly.
*
* Mathematically,
* kr = ceil(right / delta), if !isRightIncluded
* kr = floor(right / delta) + 1, if isRightIncluded
* The double division rd = right / delta never overflows and is
* exact, except in the presence of underflows. But even underflows
* do not affect the outcomes of ceil() and floor(), except,
* in turn, when the result drops to 0, that is, rd = 0.
*
* crd is a corrected version of rd when rd is zero. It is simply
* right / delta, but rounded away from 0 to preserve information
* ensuring correct outcomes in ceil() and floor().
*
* We know that -2^P <= kl, so
* -2^P <= kl + nextLong(n)
* Also, since right <= -left, we know that
* kr <= -kl + 1
* so that
* 0 < n <= -2 kl + 1
* This implies
* kl + nextLong(n) <= kl + (-2 kl) = -kl <= 2^P
* and thus
* -2^P <= kl + nextLong(n) <= 2^P
* which shows that kl + nextLong(n) can be cast exactly to double.
*
* Further, if isLeftIncluded then left = kl delta, so that we get
* left = kl * delta <= (kl + nextLong(n)) * delta
* For any other k < kl, when nextLong(n) = 0 we would have
* (k + nextLong(n)) * delta < left
* Otherwise, left = (kl - 1) delta, and therefore
* left = (kl - 1) * delta < (kl + nextLong(n)) * delta
* For any other k < kl, when nextLong(n) = 0 we would get
* (k + nextLong(n)) * delta <= left
* Either way, the lhs expression would not belong to I.
* That is, kl is the smallest integer such that kl delta always
* lies in I (it is an exact double).
*
* Similar considerations show that kr is the smallest integer such
* that kr delta lies to the right of I (it is an exact double).
*
* All the above means that (kl + nextLong(n)) * delta is an exact
* double lying in I and that kl and kr, thus n, are the best
* possible choices to ensure the largest number of equidistant
* doubles in I. Uniform distribution relies on the guarantee
* afforded by nextLong().
*/
delta = nextUp(left) - left;
double rd = right / delta;
double crd = rd != 0 || right == 0 ? rd : copySign(Double.MIN_VALUE, right);
kr = isRightIncluded ? (long) floor(crd) + 1 : (long) ceil(crd);
kl = (long) (left / delta) + (isLeftIncluded ? 0 : 1);
} else {
/* Here,
* right > 0, -right < left <= right
*
* Considerations similar to the ones above apply here as well.
*/
delta = right - nextDown(right);
double ld = left / delta;
double cld = ld != 0 || left == 0 ? ld : copySign(Double.MIN_VALUE, left);
kl = isLeftIncluded ? (long) ceil(cld) : (long) floor(cld) + 1;
kr = (long) (right / delta) + (isRightIncluded ? 1 : 0);
}
n = kr - kl;
return DoubleStream.generate(() -> (kl + nextLong(n)) * delta).sequential();
}
/**
* Returns an effectively unlimited stream of pseudorandomly chosen
* {@code int} values.

View File

@ -0,0 +1,349 @@
/*
* Copyright (c) 2023, Oracle and/or its affiliates. All rights reserved.
* DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
*
* This code is free software; you can redistribute it and/or modify it
* under the terms of the GNU General Public License version 2 only, as
* published by the Free Software Foundation.
*
* This code is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
* version 2 for more details (a copy is included in the LICENSE file that
* accompanied this code).
*
* You should have received a copy of the GNU General Public License version
* 2 along with this work; if not, write to the Free Software Foundation,
* Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
*
* Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
* or visit www.oracle.com if you need additional information or have any
* questions.
*/
import jdk.test.lib.RandomFactory;
import org.junit.jupiter.params.ParameterizedTest;
import org.junit.jupiter.params.provider.Arguments;
import org.junit.jupiter.params.provider.MethodSource;
import java.util.Iterator;
import java.util.TreeSet;
import java.util.random.RandomGenerator;
import java.util.stream.DoubleStream;
import static org.junit.jupiter.api.Assertions.*;
import static org.junit.jupiter.params.provider.Arguments.arguments;
/**
* @test
* @bug 8302987
* @key randomness
*
* @summary Check consistency of RandomGenerator::equiDoubles
* @library /test/lib
* @run junit EquiDoublesTest
*
*/
public class EquiDoublesTest {
private static final int SAMPLES = 100_000;
/*
* A factor to use in the tight*() tests to make sure that
* all equidistant doubles are generated.
*/
private static final long SAFETY_FACTOR = 100L;
private static final RandomGenerator rnd = RandomFactory.getRandom();
private static double nextUp(double d, int steps) {
for (int i = 0; i < steps; ++i) {
d = Math.nextUp(d);
}
return d;
}
private static double nextDown(double d, int steps) {
for (int i = 0; i < steps; ++i) {
d = Math.nextDown(d);
}
return d;
}
static Arguments[] equi() {
return new Arguments[] {
arguments(0.0, 1e-9),
arguments(1.0, 1.1),
arguments(1.0e23, 1.1e23),
arguments(1.0e300, 1.1e300),
arguments(-1.2, 1.1),
arguments(-1.2e-30, 1.1e6),
arguments(-Double.MIN_VALUE, Double.MIN_VALUE),
arguments(-Double.MAX_VALUE, Double.MAX_VALUE),
};
}
@ParameterizedTest
@MethodSource
void equi(double l, double r) {
double[] minmax = new double[2];
resetMinmax(minmax);
DoubleStream equi = rnd.equiDoubles(l, r, true, false);
equi.limit(SAMPLES).forEach(d -> updateMinmax(minmax, d));
assertTrue(l <= minmax[0]);
assertTrue(minmax[1] < r);
resetMinmax(minmax);
equi = rnd.equiDoubles(l, r, true, true);
equi.limit(SAMPLES).forEach(d -> updateMinmax(minmax, d));
assertTrue(l <= minmax[0]);
assertTrue(minmax[1] <= r);
resetMinmax(minmax);
equi = rnd.equiDoubles(l, r, false, true);
equi.limit(SAMPLES).forEach(d -> updateMinmax(minmax, d));
assertTrue(l < minmax[0]);
assertTrue(minmax[1] <= r);
resetMinmax(minmax);
equi = rnd.equiDoubles(l, r, false, false);
equi.limit(SAMPLES).forEach(d -> updateMinmax(minmax, d));
assertTrue(l < minmax[0]);
assertTrue(minmax[1] < r);
/* with negated intervals */
resetMinmax(minmax);
equi = rnd.equiDoubles(-r, -l, true, false);
equi.limit(SAMPLES).forEach(d -> updateMinmax(minmax, d));
assertTrue(-r <= minmax[0]);
assertTrue(minmax[1] < -l);
resetMinmax(minmax);
equi = rnd.equiDoubles(-r, -l, true, true);
equi.limit(SAMPLES).forEach(d -> updateMinmax(minmax, d));
assertTrue(-r <= minmax[0]);
assertTrue(minmax[1] <= -l);
resetMinmax(minmax);
equi = rnd.equiDoubles(-r, -l, false, true);
equi.limit(SAMPLES).forEach(d -> updateMinmax(minmax, d));
assertTrue(-r < minmax[0]);
assertTrue(minmax[1] <= -l);
resetMinmax(minmax);
equi = rnd.equiDoubles(-r, -l, false, false);
equi.limit(SAMPLES).forEach(d -> updateMinmax(minmax, d));
assertTrue(-r < minmax[0]);
assertTrue(minmax[1] < -l);
}
private void resetMinmax(double[] minmax) {
minmax[0] = Double.POSITIVE_INFINITY;
minmax[1] = Double.NEGATIVE_INFINITY;
}
private void updateMinmax(double[] minmax, double d) {
if (d < minmax[0]) {
minmax[0] = d;
}
if (d > minmax[1]) {
minmax[1] = d;
}
}
static Arguments[] tight() {
return new Arguments[] {
arguments(0.0, (short) 100),
arguments(1.0, (short) 100),
arguments(1.1, (short) 100),
arguments(1.0e23, (short) 100),
arguments(1.0e300, (short) 100),
arguments(-1.2, (short) 100),
arguments(-1.2e-30, (short) 100),
arguments(-Double.MIN_VALUE, (short) 100),
arguments(-Double.MIN_VALUE, (short) 2),
arguments(-Double.MAX_VALUE, (short) 2),
};
}
/*
* All equidistant doubles in a tight range are expected to be generated.
* The arguments must be chosen as to not overlap a value with irregular
* spacing around it.
*/
@ParameterizedTest
@MethodSource
void tight(double l, short steps) {
double r = nextUp(l, steps);
TreeSet<Double> set = new TreeSet<>();
DoubleStream equi = rnd.equiDoubles(l, r, true, false);
equi.limit(SAFETY_FACTOR * steps).forEach(set::add);
assertEquals(steps, set.size());
checkEquidistance(set);
set.clear();
equi = rnd.equiDoubles(l, r, true, true);
equi.limit(SAFETY_FACTOR * steps).forEach(set::add);
assertEquals(steps + 1, set.size());
checkEquidistance(set);
set.clear();
equi = rnd.equiDoubles(l, r, false, true);
equi.limit(SAFETY_FACTOR * steps).forEach(set::add);
assertEquals(steps, set.size());
checkEquidistance(set);
set.clear();
equi = rnd.equiDoubles(l, r, false, false);
equi.limit(SAFETY_FACTOR * steps).forEach(set::add);
assertEquals(steps - 1, set.size());
checkEquidistance(set);
/* with negated intervals */
set.clear();
equi = rnd.equiDoubles(-r, -l, true, true);
equi.limit(SAFETY_FACTOR * steps).forEach(set::add);
assertEquals(steps + 1, set.size());
checkEquidistance(set);
set.clear();
equi = rnd.equiDoubles(-r, -l, true, true);
equi.limit(SAFETY_FACTOR * steps).forEach(set::add);
assertEquals(steps + 1, set.size());
checkEquidistance(set);
set.clear();
equi = rnd.equiDoubles(-r, -l, false, true);
equi.limit(SAFETY_FACTOR * steps).forEach(set::add);
assertEquals(steps, set.size());
checkEquidistance(set);
set.clear();
equi = rnd.equiDoubles(-r, -l, false, false);
equi.limit(SAFETY_FACTOR * steps).forEach(set::add);
assertEquals(steps - 1, set.size());
checkEquidistance(set);
}
static Arguments[] tightWithIrregularSpacing() {
return new Arguments[] {
arguments(0x1p-1, (short) 15, (short) 23),
arguments(0x1p0, (short) 17, (short) 5),
arguments(0x1p1, (short) 7, (short) 8),
arguments(0x1p-600, (short) 28, (short) 33),
arguments(0x1p600, (short) 9, (short) 19),
};
}
/*
* m must be a power of 2 greater than Double.MIN_NORMAL
*/
@ParameterizedTest
@MethodSource
void tightWithIrregularSpacing(double m, short lSteps, short rSteps) {
double l = nextDown(m, 2 * lSteps);
double r = nextUp(m, rSteps);
int steps = lSteps + rSteps;
TreeSet<Double> set = new TreeSet<>();
DoubleStream equi = rnd.equiDoubles(l, r, true, false);
equi.limit(SAFETY_FACTOR * steps).forEach(set::add);
assertEquals(steps, set.size());
checkEquidistance(set);
set.clear();
equi = rnd.equiDoubles(l, r, true, true);
equi.limit(SAFETY_FACTOR * steps).forEach(set::add);
assertEquals(steps + 1, set.size());
checkEquidistance(set);
set.clear();
equi = rnd.equiDoubles(l, r, false, true);
equi.limit(SAFETY_FACTOR * steps).forEach(set::add);
assertEquals(steps, set.size());
checkEquidistance(set);
set.clear();
equi = rnd.equiDoubles(l, r, false, false);
equi.limit(SAFETY_FACTOR * steps).forEach(set::add);
assertEquals(steps - 1, set.size());
checkEquidistance(set);
/* with negated intervals */
set.clear();
equi = rnd.equiDoubles(-r, -l, true, true);
equi.limit(SAFETY_FACTOR * steps).forEach(set::add);
assertEquals(steps + 1, set.size());
checkEquidistance(set);
set.clear();
equi = rnd.equiDoubles(-r, -l, true, true);
equi.limit(SAFETY_FACTOR * steps).forEach(set::add);
assertEquals(steps + 1, set.size());
checkEquidistance(set);
set.clear();
equi = rnd.equiDoubles(-r, -l, false, true);
equi.limit(SAFETY_FACTOR * steps).forEach(set::add);
assertEquals(steps, set.size());
checkEquidistance(set);
set.clear();
equi = rnd.equiDoubles(-r, -l, false, false);
equi.limit(SAFETY_FACTOR * steps).forEach(set::add);
assertEquals(steps - 1, set.size());
checkEquidistance(set);
}
private void checkEquidistance(TreeSet<Double> set) {
if (set.size() < 3) {
return;
}
Iterator<Double> iter = set.iterator();
double prev = iter.next();
double curr = iter.next();
double delta = curr - prev;
while (iter.hasNext()) {
prev = curr;
curr = iter.next();
assertEquals(delta, curr - prev);
}
}
static Arguments[] empty() {
return new Arguments[] {
arguments(1.0),
arguments(-1.0),
arguments(0.0),
arguments(nextDown(Double.MAX_VALUE, 1)),
arguments(nextUp(-Double.MAX_VALUE, 1)),
};
}
@ParameterizedTest
@MethodSource
void empty(double l) {
assertThrows(IllegalArgumentException.class,
() -> rnd.equiDoubles(l, l, true, false)
);
assertThrows(IllegalArgumentException.class,
() -> rnd.equiDoubles(l, nextUp(l, 1), false, false)
);
assertThrows(IllegalArgumentException.class,
() -> rnd.equiDoubles(nextDown(l, 1), l, false, false)
);
assertThrows(IllegalArgumentException.class,
() -> rnd.equiDoubles(l, l, false, true)
);
assertThrows(IllegalArgumentException.class,
() -> rnd.equiDoubles(l, nextDown(l, 1), true, true)
);
assertThrows(IllegalArgumentException.class,
() -> rnd.equiDoubles(nextUp(l, 1), l, true, true)
);
}
}