Perfect Random Floating-Point Numbers
56 comments
·May 4, 2025camel-cdr
lblume
How relevant are subnormals in this context?
i.e. what percentage of the range do they occupy and how often would one expect to need to sample to obtain one?
kbolino
There's as many subnormals as there are distinct values with the same exponent, so 2^23 for floats and 2^52 for doubles. This includes the number 0, but you may consider it a special case different from the other subnormals; if so, subtract one. How much of the range they occupy also depends on the interval you're requesting. However, for the common interval from 0 to 1, they occupy about 1/127th of the range for floats and 1/1023rd of the range for doubles.
One can reason to this conclusion by noting that the subnormals have the same effective exponent as the smallest normal number (except that they use the 0.xxx representation instead of 1.xxx; the exponent in question is -126 for floats and -1022 for doubles), and then observing that all normal numbers with exponents from that minimum up to and including -1 are also in the range. Since the number 1 is usually not included in the interval, we can ignore it, but even if it was, it would be the only representative from its exponent (0) anyway, so the effect on the size of the range would be miniscule.
Dylan16807
You will never ever ever obtain one.
Even calculating your exponent with a single count_trailing_zeros on a single 64 bit value will give you results indistinguishable from perfect.
FabHK
There's another approach for doing this: generate a random number between 1 and 2 (they all have the same exponent) and then subtract 1 (that's the default implementation in Julia [0]). But I think it also just gives you 2^53 different numbers.
So, between .5 and 1 you miss out on every second representable number, between 0.25 and .5 you miss out on 3/4 of them, and so on.
I guess for many cases that's good enough, but the article seems like a nice improvement.
ETA: Lemire has some thoughts on this [1] and links to what might be a prior solution [2]. Vigna (of xoroshiro fame) writes about it at the bottom of [3] and also links to [2]. So, presumably the implementation described in the article is faster? ("There have been some past attempts to fix these flaws, but none that avoid a huge performance penalty while doing so."
EDIT2: BTW, one of the things I love about HN (well, the world, really) is that there are people that care deeply that we can uniformly sample floats between 0 and 1 correctly, and all of them, and do it faster.
[0] see https://github.com/JuliaLang/julia/blob/master/stdlib/Random...
rand(r::AbstractRNG, ::SamplerTrivial{CloseOpen01_64})
= rand(r, CloseOpen12()) - 1.0
[2] https://mumble.net/~campbell/2014/04/28/uniform-random-float
pclmulqdq
That is essentially just the fixed point method mentioned here, and it is uniformly distributed across the range (and faster because you do less), but it does not fill the low bits of the number correctly (53 bits of precision != 2^53 of them between 0 and 1). If you look at the final distribution of numbers and you use under 1-10 million numbers, that method and the division method will appear to be unbiased unless you have a pathological case.
badmintonbaseba
That looks kinda the same end result as the fixed point method, but possibly cheaper to execute.
thaumasiotes
Well, it has to be. The number can only have 53 bits of precision.
There are more than 2^53 floating point values between 0 and 1, but if you use them all, you'll fail to have a uniform distribution. Either certain regions of the space will be far more likely than other regions, or certain values will be far more likely than other values. The article takes the view that the second of those options is desirable, but it's not clear why.
And the algorithm the article describes is:
1. Generate a random float the uniform way.
2. If it's close to zero and has extra precision available, generate some more random bits to fill in the extra precision.
It's always going to be faster to skip step 2. What are we gaining by not skipping it?
Sharlin
Floats are typically used to approximate reals, a uniform distribution over individual float values is usually not what is needed. What is wanted is uniform distribution over equal-sized intervals.
2^53 equally-spaced values within [0, 1) is plenty enough for many use cases, but it's still fundamentally a set of fixed-point rather than floating-point numbers. For a true random floating-point value, all the available precision should be used such that every single value in the domain has some probability mass (except maybe subnormals, but you can usually forget about subnormals). Especially with 32-bit floats it's not that difficult to run out of precision when you start doing math and only have a fixed-point subset of size 2^23 available.
mlyle
> It's always going to be faster to skip step 2. What are we gaining by not skipping it?
Having any number be possible is nice. Having the entire space be reachable when you're sampling a solution space and the solution may be smaller than 1/8. Or having small vector quantities not have their angles ridiculously constrained by quantization.
Most applications don't need it-- indeed, a double has a lot of precision compared to what most people care about, and throwing away some of it at the low end of the randomness range is no big deal. But sometimes, a double is barely (or not even) enough.
FabHK
> certain values will be far more likely than other values
That's fine, as for the normal understanding of uniformity on the interval, you want the probability that a value is chosen to be proportional to the "space" between its left and right neighbor, and that space is not constant on the interval.
What you try to emulate is picking a real number uniformly between 0 and 1, then returning the closest floating point number (except 1, maybe).
camel-cdr
> but if you use them all, you'll fail to have a uniform distribution
Only if you generate them all with equal probability.
david-gpu
At some point in my career I worked on floating point hardware design, and I agree with your approach.
pixelpoet
Some good references on this which IMO should have been mentioned in the article:
https://pharr.org/matt/blog/2022/03/05/sampling-fp-unit-inte...
https://marc-b-reynolds.github.io/distribution/2017/01/17/De...
possiblywrong
> Second, the least significant bits that come from this generator are biased.
I don't think this is true; I'm guessing that the author borrowed this observation from one of the various other articles linked in this thread, that address the situation where we draw a 64-bit random unsigned integer and divide by `1<<64`, instead of drawing 53 bits and dividing by `1<<53`. The former does introduce a bias, but the latter doesn't (but does still limit the fraction of representable values).
camel-cdr
I'm not sure which PRNG is used here, but some PRNGs have regularities in the lower bits.
This is e.g. the case for classical LCGs and the xoroshiro*+ family of PRNGs:
> If, however, one has to generate only 64-bit floating-point numbers (by extracting the upper 53 bits) xoshiro256+ is a slightly (≈15%) faster generator with analogous statistical properties. For general usage, one has to consider that its lowest bits have low linear complexity and will fail linearity tests (https://prng.di.unimi.it/)
> When m = 2k there is a further problem. The period of the bth bit of an LCG (where bits are numbered from the right, starting at one) is 2b , thus although the period is 2k, only the high bits are good and the lower order bits exhibit a clear repeating pattern. (https://www.pcg-random.org/pdf/hmc-cs-2014-0905.pdf)
tialaramex
> For unbiased random floating-point numbers, generate floating-point numbers with probabilities given by drawing a real number and then rounding to floating point.
Immediately there are alarms going off in my mind. Your machine definitely can't pick real numbers, Almost All of them are non-computable. So you definitely can't be doing that, which means you should write down what you've decided to actually do because it's not that.
What you actually want isn't the reals at all, you want a binary fraction (since all your f64s are in fact binary fractions) between 0 and 1, and so you just need to take random bits until you find a one bit (the top bit of your fraction), counting all the zeroes to decide the exponent, then you need 52 more bits for your mantissa.
I'm sure there's a faster way to get the same results, but unlike the article's imagined "drawing a real number" this is actually something we can realize, since it doesn't involve non-computable numbers.
Edited: You need 52 more bits, earlier this comment said 51 but the one bit you've already seen is implied in the floating point type, so we need 52 more, any or all of which might be zero.
yorwba
The part you're quoting is the spec, not the implementation. It implicitly assumes a machine model where sampling a uniformly random real number in an interval and correctly rounding it to a floating-point number are both supported operations. Of course this cannot be a real machine or even a Turing machine, but you can imagine and mathematically model a more powerful machine which can do these things. (Turing machines are imaginary anyways.)
This unimplementable algorithm for an imaginary machine is still useful, because it nails down the probability distribution we expect the output samples to have. And because they're sampled from a finite value space and the probabilities of obtaining each possible output are rational numbers, we have a chance of finding an implementable algorithm for a real machine that produces the same output distribution.
The rest of the article goes into detail on what they decided to actually do. It's similar to your proposal but has more edge-case handling and deals with different rounding modes.
lblume
> Turing machines are imaginary anyways.
Sure, but there is a simple correspondence between modern/typical computers (that does require a potentially unlimited supply of new hard drives, but that is rarely relevant) and Turing machines, while such a correspondence can probably never exist for the "real number" model.
badmintonbaseba
Generalizing this to arbitrary intervals, not just [0, 1) looks tricky. Just scaling and shifting a perfect uniform result from [0, 1) doesn't result in a perfect random floating-point number.
TimWolla
Indeed. It might even result in out of range values due to implicit rounding (`3.5 + (4.5 - 3.5) * (the largest float smaller than 1)` results in 4.5, which is unexpected).
I'm a contributor to the PHP programming language and a maintainer of the randomness functionality and did some research into this topic to provide a floating point generation method that is as good as possible and came across the "Drawing Random Floating-point Numbers from an Interval" paper by Frédéric Goualard.
PHP 8.3 ships with a method generating random floats from arbitrary intervals based on the γ-section algorithm described in that paper. The PHP documentation goes into additional detail:
https://www.php.net/manual/en/random-randomizer.getfloat.php
As part of the implementation I also reached out to Prof. Goualard with some questions for clarification regarding the behavior for extreme values, which resulted in an Corrigendum being published, since the algorithm was buggy in those cases (https://frederic.goualard.net/publications/corrigendum.pdf). The entire exchange with Prof. Goualard was really pleasant.
kevmo314
Cool observation! Despite knowing both about how floating points work and how random number generation works, this never occurred to me.
I do wish the solution were a bit simpler though. Could this be upstreamed into the language's API? https://cs.opensource.google/go/go/+/refs/tags/go1.24.3:src/...
lutusp
> Could this be upstreamed into the language's API?
If a language is in use, and if people have written code that generates pseudorandom sequences based on an initial seed, then no -- bad idea. Some programs rely on the same pseudorandom sequence for a given seed, and may fail otherwise.
wtallis
That really depends on whether the language's standard library API specifies that implementations will use a particular random number generator, or merely specifies that the RNG will have certain properties. If you're going to depend on reproducible, portable RNG behavior you need to get those random numbers from an API that guarantees a particular algorithm.
hedora
I thought this was going to go further down the path that the Die Hard Battery of random number tests started:
https://www.jstatsoft.org/index.php/jss/article/view/v007i03...
(Dieharder is probably a better suite to use at this point, but that paper is approachable.)
FabHK
Better test suites are PractRand and TestU01.
FabHK
So, fun fact:
Between 0 and 1, you have about a quarter of all floating point numbers (and then a quarter > 1, a quarter < -1, and a quarter between -1 and 0).
Between 1 and 2, you have about 0.024% of all floating point numbers (for double precision, a factor of around 1023).
gitroom
Yeah I've messed with this and always end up wishing the simple ways were actually good enough. The bit where floats aren't really evenly spaced just gets annoying. Still, I kinda get the itch to cover all cases. Respect.
sfpotter
Probably relevant: https://prng.di.unimi.it/
I use the PRNGs here for my own needs and they work great... at least as far as I can tell. :-)
mahemm
Why not just read 64 bits off /dev/urandom and be done with it? All this additional complexity doesn't actually buy any "extra" randomness over this approach, and I'm skeptical that it improves speed either.
mlyle
The problem is, there's around 2^62 double precision numbers between 0 and 1, but they're not uniformly spaced-- there's many, many more between 0 and 0.5 (nearly 2^62) than there are between 0.5 and 1 (around 2^52), for instance.
So, if you want a uniform variate, but you want every number in the range to be possible to generate, it's tricky. Each individual small number needs to be much less likely than each individual large number.
If you just draw from the 2^62 space randomly, you almost certainly get a very small number.
spyrja
Seems to me that the simplest solution would be to repeatedly divide the range of numbers into two halves, then randomly selecting either one until it converges onto a single value. In C this might look something like this:
double random_real(double low, double high, int (*random_bit)(void)) {
if (high < low)
return random_real(high, low, random_bit);
double halfway, previous = low;
while (true) {
halfway = low + (high - low) / 2;
if (halfway == previous)
break;
if (random_bit() & 1)
low = halfway;
else
high = halfway;
previous = halfway;
}
return halfway;
}
seanhunter
So you have two doubles halfway and previous and a recursion that depends on if(halfway==previous) to terminate, where halfway is the result of a floating point calculation. You sure that’s going to work? I suspect it will frequently fail to terminate when you think.
Secondly, why does this generate a uniform random number? It’s not clear to me at all. It seems it would suffer the exact problem GP’s talking about here, that certain ranges of numbers would have a much higher probability than others on a weighted basis.
camel-cdr
That wouldn't produce uniformly distributed random floating-point numbers.
I've written an algorithm that generates a uniform random float in the range [a,b] that can produce every representable floating-point value with a probability proportional to the covered real number range*: https://github.com/camel-cdr/cauldron/blob/main/cauldron/ran...
* Well, almost. I messed up the probability around subnormals slightly and haven't gotten around to fixing it yet.
Still, the algorithm it self should works for all other values and it was measurably more accurate than the regular algorithm.
In terms of performance it is 10x slower compared to the regular implementation of `(randu32(x)>>8)*(1.0f/(1<<24))` on my Zen1 desktop.