This crate implements several techniques for converting streams of random bits into floating-point numbers distributed as a uniform real in the unit interval.
The division technique of extracting 53 bits and dividing by 2⁵³, used, for
example, by the rand crate, provides 2⁵³
equispaced floating-point numbers, but cannot generate almost all the
floating-point numbers in the unit interval. The techniques in this crate will
generate correctly all floating-point numbers in the unit interval, including
subnormals, with a probability proportional to the space around them.
The uniform module exposes the technique of choice (an alias for the
pekkizen module described below). Its unif_01 function turns any
source of random 64-bit words into a uniform f64 in [0 . . 1); with the
rand feature, enabled by default, the Unif01Ext extension trait
provides the same conversion as a unif_01 method on every
generator of the rand crate:
use rand_float::uniform::Unif01Ext;
// With any generator of the rand crate:
let x = rand::rng().unif_01();
assert!((0.0..1.0).contains(&x));
// With any source of random 64-bit words:
let mut src = rand_float::sources::Weyl(42);
let y = rand_float::uniform::unif_01(|| src.next_u64());
assert!((0.0..1.0).contains(&y));A relatively innocuous pull request ended up in a rabbit hole that led to this crate. The motivation was the state of disarray of the literature on the topic: academic papers, GitHub repositories, and free-floating sources, often lacking comparison with other approaches, a few bugs, and in some cases definitely hostile notation. By gathering all techniques in the same place, we have a common ground for comparison and benchmarking, and cross references for future implementations. A detailed discussion of the issues with the current techniques can be found in this paper by Frédéric Goualard.
The provided gaussian_tail example shows the statistical defects caused by the
division technique when generating a large number of Gaussian deviates by
inversion:
1.2e12 draws per converter on 10 threads, seed 42
x/2^53 + zero guard:
[...]
286733 deviates beyond -5.035σ in 118 s; 17 duplicates
-5.7127σ drawn 2 times
-5.3617σ drawn 2 times
-5.3529σ drawn 2 times
...
unif_01:
[...]
286748 deviates beyond -5.035σ in 246 s; 0 duplicates
The detected collisions are statistically impossible: their presence is only due to the fact that the division technique quantizes the tail of the Gaussian distribution, as below 2⁻²² it can return just 2³¹ distinct values, all multiples of 2⁻⁵³. For the same reason it cuts off the tail arbitrarily, as it never returns nonzero numbers closer to zero than 2⁻⁵³. Details about this specific issue can be found in this ACM Computing Survey (section 3.1).
This is also the same problem that Cleve Moler addressed in MatLab using a custom generation scheme.
Another area of application of these techniques is that of secure computations, as the lower bits of the numbers generated by the division technique are not random at all. In particular, differential privacy needs floating-point numbers that are random in every bit.
Every technique is implemented as a pure transformation of a source of
uniform random 64-bit words (any FnMut() -> u64), documented with its
exact output distribution, and benchmarked against the others.
The translation from the sources in different languages was operated by Anthropic's Fable and extensively checked against the original implementation. In the process, we isolated a few bugs in the original code, that have been reported to the authors.
| Module | Origin | Distribution | Reachable values | Words per f64 |
|---|---|---|---|---|
division |
folklore | equispaced | the 2⁵³ multiples of 2⁻⁵³ in [0 . . 1) | 1 |
pekkizen |
Pekka Pulkkinen's uniFloats (Float64_64) |
uniform real rounded down to a 2⁻⁶⁴ grid | every float in [2⁻¹² . . 1); 2⁵² values spaced 2⁻⁶⁴ below 2⁻¹² | 1 |
pekkizen |
Pekka Pulkkinen's uniFloats (Float64_117) |
uniform real rounded down to a 2⁻¹¹⁷ grid | every float in [2⁻⁶⁵ . . 1); multiples of 2⁻¹¹⁷ below 2⁻⁶⁵ | 1 + ≈2⁻¹² |
pekkizen |
Pekka Pulkkinen's uniFloats (Float64full) |
uniform real in [0 . . 1) rounded down | every float in [0 . . 1), including all subnormals | 1 + ≈2⁻¹² |
campbell |
Taylor R. Campbell's binary64fast.c |
uniform real in [0 . . 1] rounded to nearest | every float in [2⁻¹²⁸ . . 1] | 2 (or 3, for constant time) |
campbell |
Taylor R. Campbell's random_real.c |
uniform real in [0 . . 1] rounded to nearest | every float in [2⁻¹⁰²⁴ . . 1], and 0 | ≈1.5 |
badizadegan |
fp-rand (round-down variant) | uniform real in (0 . . 1) rounded down | every float in [0 . . 1), including all subnormals | 1 + ≈2⁻¹² |
A few observations:
-
Some of the techniques have variants: for example, both Badizadegan's and Pulkkinen's can return values in [0 . . 1] by rounding to nearest, but the semi-open interval is our target as it is the right one to do inversion (e.g., if you do inversion of a discrete distribution on [0 . . 1] you'll be in trouble).
-
Campbell's constant-time code is geared towards shielding from timing side-channel attacks, which is a non-goal for us.
-
The documentation of
campbell::realdiffers from the comments in Campbell'srandom_real.c: while porting it we found a minor bug (the all-zeros cutoff fires one 64-bit word too early, so the subnormals below 2⁻¹⁰²⁴ are unreachable and 0 is returned slightly too often), which we reported to the author. The code is a faithful port; our documentation describes what the code actually does. -
Badizadegan's technique is almost optimal in terms of entropy, a property shared by none of the other techniques. To use this property, however, one should keep track of the used bits in each 64-bit output, making the stateless approach of this crate impossible.
use rand_float::{badizadegan, campbell, division, pekkizen};
struct Xoroshiro128pp([u64; 2]);
impl Xoroshiro128pp {
fn next_u64(&mut self) -> u64 {
let [s0, mut s1] = self.0;
let result = s0.wrapping_add(s1).rotate_left(17).wrapping_add(s0);
s1 ^= s0;
self.0 = [s0.rotate_left(49) ^ s1 ^ (s1 << 21), s1.rotate_left(28)];
result
}
}
let mut src = Xoroshiro128pp([0x243F6A8885A308D3, 0x13198A2E03707344]);
let a = division::f64_53bits(|| src.next_u64());
let b = pekkizen::f64_64(|| src.next_u64());
let c = campbell::fast(|| src.next_u64());
let d = badizadegan::f64_down(|| src.next_u64());cargo bench drives every technique with the same Weyl-sequence source in two
settings: one conversion per call, and filling an array of 1024 doubles per
iteration.
All builds use -C target-cpu=native (set in .cargo/config.toml), so the
benchmarks measure code generated for the machine they run on.
On CPUs with heterogeneous cores, pin the run to one core (build first so compilation stays parallel):
cargo bench --no-run
taskset -c 2 cargo bench # Linux; pick a performance coreTo run the benchmarks and turn the results into a bar chart (ns per
generated f64, single-call and array-fill bars side by side; requires
Python with matplotlib) in one line:
cargo bench 2>&1 | python3 python/plot_bench.py -o bench.pdfThe output format follows the extension (.pdf, .png, .svg, ...), and
the script also accepts a previously saved log: python3 python/plot_bench.py bench.txt -o bench.pdf.
Here we display the results of the benchmarks on some architectures. Note that because of a quirk in LLVM, we had to add (at least for the time being) a cold barrier preventing some un-optimization of the array case due to interference with the underlying Weyl generator. The barrier sometimes however disturbs a bit the single-call case. We hope to remove it in the future.
I would like to thank Nima Badizadegan, Taylor R. Campbell, Frédéric Goualard, and Reiner Pope for interesting discussions and a lot of useful pointers.
Original code in this crate is dual-licensed under Apache-2.0 or MIT. The ported techniques retain the licenses of their authors, reproduced in the header of the respective source files: