Uniformly distributed secure floating point numbers in [0,1)

Question

Is there any way to generate

• cryptographically secure
• uniformly distributed

floating point numbers in the interval [0,1)?

For example, in Javascript, there is Math.random(), which yields these (pseudo-random) numbers, yet not cryptographically secure.

• There are several algorithms at https://stackoverflow.com/questions/686353/c-random-float-number-generation, yet at https://stackoverflow.com/questions/34575635/cryptographically-secure-float, the problem of denormal numbers was pointed out.
• Inversion sampling as of Wikipedia requires a uniformly distributed [0,1) number to start, so this is seemingly not trivial to apply.

What is a good approach?

Answer 1

For most purposes, the following trivial method should be adequate:

1. Generate a random $n$-bit integer (i.e. a uniformly distributed random integer between $0$ and $2^n-1$ inclusive).

2. Cast the random integer into a float and multiply it by $2^{-n}$ to scale it to the interval $[0,1)$.

There's some leeway in choosing the number of bits $n$, but a generally reasonable choice is the significand precision $n_p$ of your floating-point format (e.g. 53 bits for IEEE 754 double-precision floats). In particular, this is the highest number of bits for which the distribution of the resulting floats will actually have the full $n$ bits of entropy — or, in other words, the highest number of bits for which step 2 above is fully reversible. If $n > n_p$, the conversion to float will map some pairs of integers to the same floats, resulting in an output distribution with less than $n$ (but still more than $n_p$) bits of entropy.

Regardless of the value of $n$ chosen, the distribution of the floats generated using this method will be "approximately uniform", in the sense that, for any $0 \le a < b < 1$, the probability of the output falling in the subinterval $[a,b)$ will differ from $b - a$ by at most $\max(2^{-n}, 2^{-n_p})$. For $n \le n_p$, it will also be uniform in the sense that each float value that can be generated by this method is generated with the same probability ($2^{-n}$).

However, the floats generated using this method will not be uniform in the sense that each possible float within $[0,1)$ would be equally likely to be generated. In fact, this holds for any random float generation method that is approximately uniform in the sense described above, simply because the density of representable floating point values is itself not uniform — there are far more floats in the interval $[0,\frac12)$ than in $[\frac12, 1)$. Thus, it's fundamentally impossible to generate random floats that would be uniformly distributed in both of these senses.

In any case, it's important to realize that a random float in $[0,1)$ is not the same thing as a uniformly distributed random real number in $[0,1)$: the latter distribution has infinite entropy, whereas the entropy of any actual distribution of floats representable on a computer is inevitably finite. While it's sometimes (I would even say most of the time, at least outside crypto) fine to substitute one for the other, it's important to keep the difference in mind when analyzing the security of algorithms that, conceptually, assume access to uniformly distributed real numbers.

Answer 2

tl;dr Since you asked about JavaScript specifically, here's the quick answer for the standard stackoverflow copypasta technique of modern coding:

'use strict'

function uniform01 () {
    function uniform32 () {
        let a = new Uint32Array(1)
        window.crypto.getRandomValues(a)
        return a[0]
    }

    // sample e from geometric(1/2) by counting zeros in uniform bits
    // but stop if we have too many zeros to be sensible
    let e = 0, x
    while ((x = uniform32()) == 0) {
        if ((e += 32) >= 1075) {
            return 0
        }
    }

    // count the remaining leading zeros in x
    e += Math.clz32(x)

    // sample s' = sl + sh*2^32 from odd integers in (2^63, 2^64)
    // (beware javascript signedness)
    let sl = (uniform32() | 0x00000001) >>> 0
    let sh = (uniform32() | 0x80000000) >>> 0

    // round s' to floating-point number
    let s = sl + sh*Math.pow(2, 32)

    // scale into [1/2, 1]
    let u = s * Math.pow(2, -64)

    // apply the exponent e
    return u * Math.pow(2, -e)
}

WARNING! This code does not run in constant time. In particular, the time it takes is proportional to the number of leading zeros in the fractional part of the significand. If you want it to run in constant time, you must adapt the loop to always run in the maximum number of iterations it might take. And you certainly can't use JavaScript, so you'll have to adapt it to something else anyway!

Why not the obvious uniform64()/2^64? You could sample a 64-bit integer in $[0, 2^{64})$ uniformly at random, and then divide by $2^{64}$, but then you exclude all floating-point numbers in $[0, 2^{-64})$, of which there are a great many, and which you should get with low but not negligible probability, $2^{-64}$. To put that magnitude into perspective, look up the Bitcoin network's hash rate, and consider the probability of winning a block at the current difficulty—which the Bitcoin network does every ten minutes—after a single try.

Background. What is the uniform distribution on floating-point numbers in $[0,1]$? The uniform distribution on real numbers in $[0, 1]$ is given by the Lebesgue measure determined by $\mu([a, b]) = b - a$. The natural definition of the uniform distribution on floating-point numbers in $[0, 1]$ for a given rounding map $\rho\colon \mathbb R \to \mathbb{FP}$ is the pushforward measure $\rho_*\mu\colon S \mapsto \mu(\rho^{-1}[S])$. Since $\mathbb{FP}$ is discrete, this has the probability mass function $x \mapsto \mu(\rho^{-1}(x))$; that is, the probability of a particular floating-point number $x$ in this distribution is the measure of real numbers in $[0,1]$ that are rounded by $\rho$ to $x$.

Technique. How do we sample from this distribution? Suppose we represent the real numbers in $[0,1]$ by their infinite binary expansions $$0.b_0 b_1 b_2 b_3 \ldots = b_0/2 + b_1/4 + b_2/8 + b_3/16 + \cdots,$$ where $0.1111\dots = 1$ as usual. For a uniform real random variable $X$ in $[0,1]$, the probability that $0 \leq X \leq 1/2$ is $1/2$, and likewise $1/2 \leq X \leq 1$, so the bit $b_0$ is a fair coin toss. Likewise, the probability that either $0 \leq X \leq 1/4$ or $1/2 \leq X \leq 3/4$ is also $1/2$, so the bit $b_1$ is also a fair coin toss, and so on. This suggests the procedure of (a) sampling an infinite sequence of bits $b_0, b_1, b_2, \dots$ uniformly at random, and then (b) rounding the real number $0.b_0 b_1 b_2 \dots$ to a floating-point number.

First approximation. Obviously we can't actually sample an infinite sequence of bits uniformly—but after a certain finite number of bits, the answer is guaranteed never to change, because there is a limited number of bits after the ‘binary point’ in any finite floating-point system. So you might be tempted to just sample an integer in $[0, 2^{1100})$ and divide by $2^{1100}$, if the prospect of reliable bignum arithmetic and floating-point conversions in cryptographic applications weren't enough to make you lose your lunch. Fortunately, it doesn't have to be that bad.

Second approximation. After the first $1$ bit, there are guaranteed to be no more $1$ bits beyond the next $p - 1$ bit positions, where $p$ is the precision of the floating-point format, i.e. the number of bits in the significand. Here $p = 53$ for IEEE 754 binary64 (‘double-precision’) floating-point. So you really need only sample bits until you get a $1$ bit, counting all the zeros toward the exponent, and then sample a handful additional bits for the significand—certainly no more than 64 additional ones, for binary64 floating-point.

Ties. For the default round-to-nearest/ties-to-even rounding map, there's a catch. The set of real numbers triggering the ‘ties-to-even’ rule has measure zero, but if we sample only finitely many bits, then the probability of triggering the rule will be nonzero, leading to a bias toward ‘even’ floating-point numbers—i.e., toward floating-point numbers whose least significant significand bit is $0$. You can trick this by always forcing the rule not to be invoked by drawing more than $p$ additional bits and setting the least significant one, which will be rounded away, to be $1$.

Summary.

1. Sample bits uniformly at random, counting $0$ bits into an exponent $e$, until you get a $1$ bit. You can stop at 1075 $0$ bits: your computer is broken and/or you've won the lottery many times over and been struck by lightning, etc. In any case, if you got that many $0$ bits, the result will be rounded to zero anyway no matter how many more coins you put in the slot machine, you poor dopamine addict.
2. Sample an additional 64 bits $s$ uniformly at random.
3. Set $s'$ to be $s$ with the high bit set to represent a leading $1.$ before the binary point in $1.b_0 b_1 b_2 \dots \times 2^{\mathrm{exponent}}$ notation, and with the low bit set to break ties, so that $s'$ is a uniform random choice of odd integer in $(2^{63}, 2^{64})$.
4. Round $s'$ to the nearest floating-point number $\hat s \in [2^{63}, 2^{64}]$. (There will never be a tie!)
5. Compute the floating-point number $u = \hat s / 2^{64} \in [1/2, 1]$. This division is exact because the significand is preserved, while the exponent of $\hat s$ was $63$ or $64$ and of $u$ is therefore $-1$ or $0$.
6. Compute and yield the floating-point quotient $\operatorname{round(u/2^e)}$. This division may not be exact, but only because the exponent may drop below the minimum exponent and be rounded to a subnormal or zero.

Make sure to use a cryptographic bit sampler, of course—use window.crypto.getRandomValues if you're confined to the dungeon of JavaScript.

$0$ and $1$. You said $[0, 1)$, but I said $[0, 1]$ throughout this answer. Why? The rounding map may send real numbers in $[0, 1)$ to the floating-point number $1$. For example, $1 - 2^{-55}$ will be rounded to $1$ in binary64 floating-point—as will any other real number more than halfway from $1 - 2^{-53}$ to $1$. The probability of rounding a uniform random choice of real number in $[0, 1]$ to $1$, namely $2^{-54}$, is small, but very much not negligible in cryptography. Do you really mean to work in $[0, 1 - 2^{54}) \subset \mathbb R$, or do you want a floating-point approximation to $[0, 1)$, which might involve the floating-point number $1$?

If you can't handle $1$, you could always do rejection sampling—or, motivated by the nonuniformity of the density of floating-point numbers in $[0,1]$, you could choose a different space in which to work, such as log-space or log-odds-space, but that is a lesson for another day.

On the other talon, in binary64 floating-point, the probability of getting the floating-point number $0$ is $2^{-1075}$ (or a whopping $2^{-1023}$ if subnormals are flushed to zero), which in cryptographic terms is negligible to the fourth power at a 256-bit security level! Even in binary32 (‘single-precision’), it's $2^{-150}$ (or $2^{-127}$ in flush-to-zero mode). So you never have to worry about getting $0$ from a serious uniform floating-point $[0, 1]$ sampler in practice.

Unfortunately, most of the alleged uniform floating-point $[0, 1]$ samplers readily at hand in common libraries like numpy, aside from using noncryptographic PRNGs, have lurking zero traps to catch unwary numerical cryptanalysts long after your code passing all tests on millions of trials has shipped into production. Caveat iactator: sampler beware!

Subnormals. Proving whether this does the right thing even with subnormals is left as an exercise for the reader. Quantifying the probability that you even have to care is left as an easier exercise for the reader. Hint: You don't.