rand-float

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.

How to use

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));

Motivation

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.

Implementations

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.

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::real differs from the comments in Campbell's random_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.

Using specific techniques

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());

Benchmarks

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 core

To 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.pdf

The 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.

Results

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.

AMD Ryzen 9 5950X

12th Gen Intel® Core™ i7-12700KF @3.60 GHz

Intel® Xeon® X5660 @2.80 GHz

Apple M1 Max @2.50 GHz

Apple M5 Max @4.00 GHz

Acknowledgments

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.

Licensing

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:

• src/badizadegan.rs — MIT, Copyright (c) 2025 Nima Badizadegan (fp-rand);
• src/campbell.rs — BSD-2-Clause, Copyright (c) 2014-2026 Taylor R. Campbell (from binary64fast.c and random_real.c);
• src/pekkizen.rs — Copyright (c) 2020 Pekka Pulkkinen, distribution permitted (uniFloats);