Higher Quality Random Floats

Much has been written about how to generate random 64-bit integers; people might argue over which particular (CS)PRNG is best suited for a given task, but there is agreement on the general shape of the solution: keep some state around, and then extract 64 bits of randomness at a time, where each bit is equally likely to be either zero or one. This gives nice uniformly-distributed integers over the range [0, 264).

There is less agreement on how to generate random floating-point values. It is tempting to aim for what looks like a simple transformation of an integer generator: uniformly-distributed floats over the range [0, 1). There are various ways of doing such a transformation; one exposition does it as:

double pcg64_random_double(pcg64_t *rng) {
return (double)(pcg64_random(rng)>> 11) * 0x1.0p-53;
}

Whereas LuaJIT does it as:

uint64_t lj_prng_u64d(PRNGState *rs) {
uint64_t z, r=0;
TW223_STEP(rs, z, r)

return (r & 0x000fffffffffffffull) | 0x3ff0000000000000ull;
}

U64double u;
double d;
u.u64=lj_prng_u64d(rs);
d=u.d – 1.0;

And a Go version by Lemire does it as:

func toFloat64(seed *uint64) float64 {
x :=splitmix64(seed)
x &=0x1fffffffffffff
return float64(x) / float64(0x1fffffffffffff)
}

The [0, 1) output range is less useful than it might initially appear. It is often stated that [0, 1) can be transformed to [A, B) by multiplying by B-A then adding A, but this is only true for some values of A and B; for other values, a number just less than 1 can be transformed and end up equal to B due to floating-point rounding, i.e. [0, 1) is transformed to [A, B]. The possibility of a zero output is also unhelpful if the result is going to be log-transformed (for example as part of a Box-Muller transform), as log(0) explodes.

An output range of [0, 1] would be better behaved under most additive and multiplicative transforms. After such a transform, a range like [A, B] can be easily turned into [A, B) via rejection sampling, should the half-open range be desired. It will also turn out to be very easy (on either theoretical or practical grounds) to exclude 0 as a possible output, giving the log-friendly output range (0, 1].

Before trying to generate a random float in the range [0, 1], we should instead consider the easier problem of generating a random float in the range [0.5, 1]. There are 252+1 different IEEE-754 doubles in this range. The “rounding basin” of some double d in that range can be defined as the range of infinite-precision real numbers in the range [0.5, 1] that would round to d (assuming round-to-nearest). Taking ε=2-54, most of these doubles have a basin of d-ε through d+ε. The exceptions are the extremes of the range; 0.5 has a basin of 0.5 through 0.5+ε, and 1 has a basin of 1-ε through 1. In other words, there are 252-1 values in the range whose basin is 2ε wide, and 2 values in the range whose basin is only ε wide. For a uniform distribution, the probability of seeing d should equal the width of the rounding basin of d. Given all this, there is a fairly simple monotonic transform from the uniform integer range [0, 253) to the uniform float range [0.5, 1]:

Integer(s)Resultant floating-point double (ε=2-54)
00.5
1, 20.5 + 2ε
3, 40.5 + 4ε
5, 60.5 + 6ε
⋮⋮
253-7, 253-61 – 6ε
253-5, 253-41 – 4ε
253-3, 253-21 – 2ε
253-11

This transform can be expressed in code like so:

double rand_between_half_and_one() {
double d;
uint64_t x=rand_u64()>> 11;
x=((x + 1)>> 1) + (1022ull> 11;
x=((x + 1)>> 1) + (1021ull> 11;
uint64_t e=1022;
do {
if (rand_u64() & 1) break;
e -=1;
} while (e> 1022-75);
x=((x + 1)>> 1) + (e>=1;
if (–nbits==0) bits=rand_u64(), nbits=64;
e -=1;
} while (e> 1022-75);
x=(((x>> 11) + 1)>> 1) + (e=0) e=__builtin_ctzll(rand_u64());
x=(((x>> 11) + 1)>> 1) – ((e – 1011ull)=0 {
e=bits.TrailingZeros64(splitmix64(seed))
}
x=(((x>> 11) + 1)>> 1) – ((uint64(int64(e)) – 1011)
>>> Read full article>>>
Copyright for syndicated content belongs to the linked Source : Hacker News – https://www.corsix.org/content/higher-quality-random-floats