Generating all floating-point numbers

[vigna]

While the division strategy used by this crate (take 53 bits, divide by 2⁵³) provides 2⁵³ equispaced points in [0 . . 1), it is unable to generate most of the possible floating-point values—in particular, values closer to zero than 1/ 2⁵³.

This is a known problem—conversion to floating point has been a sore point in all major programming languages for a long time (Go used 63 bits, and old Java versions used 54, resulting in a nonuniform distribution, for example).

While the definition of uniformity induced by the definition above is reasonable (the probability of each number is proportional to the empty space around it), it is possible to generate all floating-point numbers (including subnormals) with a probability proportional to the space around it (that should be formalized in terms of ULP, but I'm keeping it simple). For ordinary computation this is kinda hair splitting, but for some computations (e.g., generating normal deviates) having a more accurate representation around zero is important.

After some discussion with the most known experts in the field (Nima Badizadegan, Frédéric Goualard, and Taylor R. Campbell) I put together a crate comparing the three best techniques I could put my hands on.

I think rank should provide (maybe as a subcrate?) at least one method providing full floating-point generation for demanding scientific environments. I believe that the for the standard, off-the-shelf case 53-bit division is ok, but Rust is a technically very precise language, and in this moment its handling of floating-point generation could be improved.

[dhardy]

This isn't a new topic: see #531, #1287, #1346.

I'm not saying that we made the right choices, but we ended up where we are because (1) most people don't need the extra precision (which requires extra random samples in the rare case that the result is close to zero), (2) extra precision is often lost anyway by transformations away from [0..1) and (3) after much deliberation, we opted not to implement both low- and high-precision samplers (#494).

Perhaps given how many times this has come up we should reconsider.

[vigna]

I have to apologize—I actually had added an UPDATED line pointing at the previous issues, and somehow it didn't get in.

The point is not moot, tho. For instance, this ACM survey in Section 3.1 discusses exactly this issue and like the lack of precise generation around 0 affects generation of samples from the tails of a normal distribution.

Rust's rand crate seems to be aimed to provide best-of-the-breed first, with downgrade to simpler solution if time of space are an issue. Given that the stock generator is a high-quality stream cipher, one would expect that high-quality floating-point generation would be at least available, too.

It can also stay in my crate forever—there's little difference for the user. LLVM will inline everything anyway.

For me it was an interesting trip because I was able to gather all the literature in a single place, issue a number of bug reports on paper claims and in code that are being fixed, so the whole thing is a bit more stable than it was, and that was useful independently of the interest of this crate. But, just in case you change your mind, the grunt work has already been done.

[dhardy]

the stock generator is a high-quality stream cipher

The motivation for this is more to limit harm in cases where users miss-assess their security requirements / assume without checking that rand::rng() etc. are secure. We've had multiple issues like #1793 even (people assuming for some reason that SmallRng should be a CSPRNG).

We already use multiple uniform-ranged-value sampling algorithms; incorporating yet another algorithm into the API somehow would be a mess.

BTW what is the baseline in your benchmarks?

CC @vks

[vigna]

We already use multiple uniform-ranged-value sampling algorithms; incorporating yet another algorithm into the API somehow would be a mess.

Well... that depends on your notion of "uniform". The uniform distribution on a finite internal has zero mass at any point. So you're using your own definition of "uniform". I'm not saying it's bad—it's the first definition I discuss in my PRNG page. But it's not the most accurate definition, because you will never generate the vast majority of the possible floating-point numbers.

BTW what is the baseline in your benchmarks?

Division, what this crate uses now.

[vigna]

BTW, I moved the crate from a "comparison" crate to a proper library with a Unif01Ext extension trait that makes available a unif_01 method for all Rng implementations, using the fastest technique. I think that's entirely sufficient for people interested in fine-tuned generation. I'll keep all the other techniques so the crate can act as a hotspot for people looking for these kind of techniques, possibly adding and benchmarking new ones.

[vigna]

Just to make clear the applications, I added to the crate an easy example of impossible statistical behavior that one obtains from rand floating-point generation, but not from rand-float.