Part 3: https://ericlippert.com/2019/02/07/fixing-random-part-3/
uint64 x = 0x3FF0'0000'0000'0000 |
return *(double*)&x - 1.0
If you really want you could tack on:
ret += (entropy >> 63) * 0x1.0p-53; // 1.110223024625156540e-16;
ret += (entropy >> 52) * 0x1.0p-64; // 5.421010862427522170e-20;
There is uniform absolute precision on each interval among [0.5, 1), [0.25 0.5), [0.125, 0.25), &c., but as you get closer and closer to 0, there are still just as many fraction bits available (same relative precision), so these numbers are much more precise in absolute terms (i.e. relative to 1.0).
You can generate a random number in [1, 2) by just filling in the fraction part with random bits, because the exponent is always the same. Then you can subtract 1 to get a number in the range [0, 1). But for a number in e.g. the range [0.125, 0.25), the last 3 bits of the fraction will always be 0, meaning you are skipping 7/8 of the possible floating point numbers in that range.
The smallest possible denormal number among IEEE doubles is 2^(–1074), but the smallest non-zero value you will get with this method is 2^(–52).
If you want to generate any possible floating point number in the range [0,1) with probability proportional to the difference between it and the next representable number, it gets pretty tricky / computationally expensive.
Edit: I think for the general case you want something like:
uint64 x = rand64() & 0x001F'FFFF'FFFF'FFFF
double xf = .5 * 1p-52 * x
// return xf; /* 53 bits of entropy */
if(xf < .5)
y = rand64() & 0x000F'FFFF'FFFF'FFFF
double yf = .5 * 1p-52 * 1p-52 * y
xf += yf
if(xf < .5 * 1p-52) ... /* loop about 20 times */
uint64 x = rand64()&0x000F'FFFF'FFFF'FFFF
return *(double*)&x * 4.49423283715579e+307
That multiplier seems to be essentially 2^1022, in which case there may be no quantization noise possible in multiplying by it. It may simply cause the exponent to be increased and the result normalised.
The limitation is that floating point numbers have different exponent depending on where in the range [0,1) they are, while every float in [1,2) has a exponent of 3FF. Generating a number in the latter range and subtracting 1.0 lets the floating point hardware do the work of figuring out the exponent.
On futher testing, though, it looks like you can just do ( rand53() * 1p-53 ), which will be hands down better as long as your floating-point hardware has a good int64->double conversion.
There's an explanation in the comments, and a link to the original fuller explanation and source, here:
I.e. I'm pretty sure your approach gives a non uniform distribution, but happy to discuss in more detail if you want.
(EDIT: replace "omit" by "leave out", because it's too close to "emit")
Unless I made some significant error, it is a uniform distribution over multiples of 2^-52. It's also a discrete approximation of the real uniform distribution in that for any subdivision of [0,1) into N equal (half-open) intervals, the most-populated interval has at most one more inhabitant the least populated one.
It is however strictly inferior to rand53()*2^-53 (NextDoubleInner at your github link) unless your int64->double or float-multiply operations are quite bad.
See this discussion here:
The redzen library contains RNGs that address all of the known issues, e.g. see:
That Xoshiro256StarStarRandom source has also been merged into Math.Net.
/scratches head at FAIC/
This is part three of a very long series that will explore the connections between LINQ, probability distributions, and the difficulties of sampling from arbitrary distributions.