# Output distribution of rand(Float32) and rand(Float64), thread 2
**URL:** https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184
**Category:** Internals & Design
**Tags:** float, random
**Created:** [October 19, 2023, 1:10pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184 "2023-10-19T13:10:07Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)
#### Post date: [October 19, 2023, 1:10pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/1 "2023-10-19T13:10:07Z")
A recent [hackernews post](https://news.ycombinator.com/item?id=37926005) lead me to necro [https://discourse.julialang.org/t/output-distribution-of-rand/12192](https://discourse.julialang.org/t/output-distribution-of-rand/12192). If we want to reopen that can, I think a new thread is in order.
### Summary
`rand()`, and much more acutely `rand(Float32)` generate a uniformly distributed random floating point number in the unit interval. In actuality, these functions generate a uniformly distributed fixed-point number in the unit interval, and convert it to floating point. This means that the output distribution close to zero is quite counter-intuitive, and very much not what one would naively expect from “random floating point number in the unit interval”.
This is a footgun. Do we want to do something about it? What should we do about it?
### Fixed point vs floating point numbers
Fixed point numbers have the property that adjacent numbers have constant distance. Floating point numbers have the property that adjacent numbers have roughly constant quotients. In the unit interval, this means that representable floating point numbers are much denser close to zero. For floating point numbers, one typically measures errors in quotients, not differences (relative error, in order to be scale-free; if one cared about absolute errors only, there is conceptually no point in using floating point numbers in the first place!).
If one adds a small and a large floating point number, then the excess precision of the small number is lost. This is not a problem in itself – it creates neither large relative nor absolute errors.
If one subtracts two numbers of the same rough magnitude, then this is called “cancellation”, and it can be highly problematic: The result has still small absolute error, but can have very large relative error. If we happened to use floating point numbers for any other reason than convenience or CPU support, then this is called catastrophic cancellation.
### How julia generates random floats
Julia generates random floats by first generating a random number `r = rand12()` in the interval `[1, 2)`. In this interval, floating point numbers are the same as fixed point numbers (literally!). This is therefore very easy to do and fast – just fill the significand/mantissa with a random bitstring.
Next we subtract `1`, 🎉
```julia-auto
rand() = rand12() - 1
```
At this point we should be very afraid: `rand12() - 1` contains a catastrophic cancellation if our random number in `[1,2)` happened to be close to 1. In other words, we see that our random numbers are very bad for results close to zero.
### When does it matter
This does not matter at all if we add our random numbers to something nonzero of constant size.
It does matter if we e.g. compute `1/rand()`. Or to give a less contrieved example: Suppose we want to sample from a random distribution that gives numbers in `[1, Inf)`; say from a power law, with probability density `~ 1/x^a` with parameter a\>1. The way this is typically done is by considering the inverse cumulative distribution function, `F(u) = x` such that `probability(sample() > x) = u`, where `u \in [0,1]`. For the power law, this means `u ~ 1 / x^(a - 1)`, or `F(u) ~ 1 / u ^ (1/(a-1)) = pow(u, -1 / abs(a-1))`. In other words, we will have code like e.g. with a = 2:
```julia-auto
sample() = constant / rand()
```
This is code that looks perfectly fine from a precision standpoint. We have a singularity at zero; but this is not an issue for floating point numbers. If we decided to instead write the real-arithmetic equivalent
```julia-auto
sampleBad() = constant / (1 - rand())
```
i.e. decided to place the singularity of our mapping at `1` instead of `0`, then alarm bells would start ringing – this is catastrophic cancellation. The statistics of extreme events will be completely wrong. But extreme events are one of the most important things when looking at fat-tailed distributions like power laws!
And, with the current julia random number generation, this is what we do:
```julia-auto
sample() = constant / (rand12() - 1)
```
In other words, the output of `rand()` is by construction unusable for applications that care about precision close to zero.
### An extreme example: Exact zeros in Float32
Consider: Take a uniform random real number `u \in [0,1]`, and ask: What is the probability that `convert(Float32, u) == 0.0f0`?
Well, a good guess would be `nextfloat(0.0f0) == 1.0f-45 `, but depending on rounding mode intricacies something like roughly `1.0f-37` is also a valid guess (smallest non-subnormal positive float).
It does not matter. The age of the universe is of rough order `1.0f27` nanoseconds, i.e. the correct answer should be “never observed in all past or future computations done by mankind”.
Will you observe `convert(Float32, rand(Float64)) === 0.0f0`? probably not today. Will you observe `rand(Float32) === 0.0f0`?
yeah, all the time:
```julia-auto
julia> function func()
N = 0
while true
N += 1
rand(Float32) == 0 && return N
end
end
func (generic function with 1 method)
julia> func()
22098446
julia> func()
3274928
julia> func()
19735810
```
It is egregiously misleading to say that `rand(Float32)` generates a uniform distributed random floating point number from the unit interval.
### What to do about it: Nothing
We can chose to ignore the problem. We can say: This is well-known among experts and they know how to work around it; this is how everyone does it; almost nobody is complaining.
### What to do about it: Document it
One way to approach this is to document this footgun and provide an alternative function that does the right thing ™
An example is `log`, which links to `log1p` for avoiding cancellation when evaluating logs of numbers close to 1. So we could have a `rand0p` that is linked in the docs, and it is up to users to know / read the docs to use the right random function, on penalty of incorrect statistics of extreme events close to zero.
### What to do about it: Fix it for real
The complementary approach to documenting this would be to provide an alternative `fastrand` function that produces the current counter-intuitive distribution, and change `rand` to produce the Right Thing ™
Thus we would replace a correctness footgun with a performance footgun: it is up to users to know / read the docs to use the right random function, on penalty of mysterious performance losses.
This has large compat issues (is there code that relies on the counter-intuitive cancellation?) and performance issues, depending on how we implement that: how fast can we generate the Right Thing? Can we SIMD it? With which instruction sets (i.e. does it need avx512 or is there a avx2 / neon way of doing this right)? On GPU? Will this lead to problems when people implement/port algorithms from papers that assume random numbers to be fixed-point-distributed instead of floating-point-distributed?
This is a can of worm that requires lots of discussion.
### What to do about it: Document it, partially fix it
We can split it in the middle: I.e. do something reasonably fast that generates a distribution that is less egregiously misleading, and provide document and link two alternatives `fastrand` and `rand0p`. An example would be: Exclude zero from `rand(Float64)` output but otherwise keep it, and implement `rand(Float32) = convert(Float32, rand(Float64))`.
### What do other languages do?
Java generates the same flawed distribution as julia. Numpy generates the same flawed Float64 distribution as julia, cf [numpy/numpy/random/\_common.pxd at 85d956917903624289c451bd0fcfa4743ec86988 · numpy/numpy · GitHub](https://github.com/numpy/numpy/blob/85d956917903624289c451bd0fcfa4743ec86988/numpy/random/_common.pxd#L67) , I think they do something better for Float32, though?
Shockingly zig does the right thing (or something close to it), cf [Zig](https://github.com/ziglang/zig/pull/10428/). In fact it was this link on the related hackernews discussion that renewed my hope that something can be done at all!
---
### Author:  [@adienes](https://discourse.julialang.org/u/adienes)
#### Post date: [October 19, 2023, 1:25pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/2 "2023-10-19T13:25:15Z")
what does Rust do?
---
### Author:  [@Syx\_Pek](https://discourse.julialang.org/u/Syx_Pek)
#### Post date: [October 19, 2023, 1:38pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/3 "2023-10-19T13:38:58Z")
I don’t see why fixing it is not possible and done by an external specialised package for this purpose, rather than changing how `rand()` behaves.
---
### Author:  [@adienes](https://discourse.julialang.org/u/adienes)
#### Post date: [October 19, 2023, 1:39pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/4 "2023-10-19T13:39:36Z")
I suspect “possible” in this case refers more to social obstacles than technical
---
### Author:  [@rfourquet](https://discourse.julialang.org/u/rfourquet)
#### Post date: [October 19, 2023, 1:39pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/5 "2023-10-19T13:39:51Z")
Actually, specifically for `Xoshiro`, we generate floats slightly differently ([https://github.com/JuliaLang/julia/blob/4ef353c6e4ad5156ad65cb40c3a2a61fbd324b59/stdlib/Random/src/Xoshiro.jl#L296](https://github.com/JuliaLang/julia/blob/4ef353c6e4ad5156ad65cb40c3a2a61fbd324b59/stdlib/Random/src/Xoshiro.jl#L296)), by generating uniformly an integer in `[0, 2^53-1]`, and multiplying it by `2.0^-53`. How does this approach compare to “get random float in `[1, 2)` and subtract 1”, or to the slower but more correct approach?
---
### Author:  [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)
#### Post date: [October 19, 2023, 1:49pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/6 "2023-10-19T13:49:05Z")
> [@Syx\_Pek](#):
>
> I don’t see why fixing it is not possible and done by an external specialised package for this purpose, rather than changing how `rand()` behaves.
The issue is the hidden footgun. If you are a person who would never write `1/rand()` or `log(rand())` or let it pass in code review, because you know in your heart that it suffers from catastrophic cancellation, and if all your upstream (i.e. authors of packages you use) are the same kind of expert, then there is no issue for you and there is nothing to fix.
> [@rfourquet](#):
>
> by generating uniformly an integer in `[0, 2^53-1]`, and multiplying it by `2.0^-53`. How does this approach compare to “get random float in `[1, 2)` and subtract 1”
It is fundamentally the same approach: It samples uniformly from a set of equispaced candidates. Floating point numbers are not equispaced, they have a lot of room close to zero, that’s the point of floating the point (sorry for the pun).
---
### Author:  [@mbauman](https://discourse.julialang.org/u/mbauman)
#### Post date: [October 19, 2023, 1:50pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/7 "2023-10-19T13:50:16Z")
That’s roughly equivalent, as I understand it, as it returns 2^{53} evenly spaced unique values between 0 and 1. The thing that makes the approach in [Output distribution of rand()](https://discourse.julialang.org/t/output-distribution-of-rand/12192) (thread 1) and the linked blog post clever is that it returns 2^{53} evenly spaced unique values in [0.5, 1), 2^{53} twice-as-closely-spaced unique values in [0.25, 0.5) (with half the probability), and so on…
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [October 19, 2023, 1:50pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/8 "2023-10-19T13:50:33Z")
they aren’t quite the same. the multiply approach gets 1 extra bit of randomness
---
### Author:  [@CameronBieganek](https://discourse.julialang.org/u/CameronBieganek)
#### Post date: [October 19, 2023, 2:55pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/9 "2023-10-19T14:55:41Z")
> [@foobar\_lv2](#):
>
> ### What to do about it: Nothing
>
> We can chose to ignore the problem. We can say: This is well-known among experts and they know how to work around it; this is how everyone does it; almost nobody is complaining.
Newbie question:
Does Distributions.jl make the proper corrections for the `rand()` behavior when sampling from, e.g. an exponential distribution? (Or for whichever distributions it is relevant.)
---
### Author:  [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)
#### Post date: [October 19, 2023, 3:08pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/10 "2023-10-19T15:08:46Z")
> [@adienes](#):
>
> what does Rust do?
As far as I can tell, they have the same issue, cf [https://github.com/rust-random/rand/blob/9a02c819cc1e4ec6959ae25eafbb5cf6acb68234/src/distributions/float.rs](https://github.com/rust-random/rand/blob/9a02c819cc1e4ec6959ae25eafbb5cf6acb68234/src/distributions/float.rs)
Rust gets a big 👍 because the source file documents quite clearly what they are doing. I know no rust and had to read no code, the comment explains it all.
They get a small 👎 because Open01 is even more egregiously misleading for `f32` if you don’t read the fine print: Using Open01 implies that you care about not getting an exact zero result, presumably because you want to apply a mapping with singularity at zero. Well, you’re not getting a zero, so you won’t get a warning shot with crashes / Inf / NaN appearing in your statistics. But the statistics of extreme events are still egregiously counter-intuitive, just silently.
> [@CameronBieganek](#):
>
> Does Distributions.jl make the proper corrections for the `rand()` behavior when sampling from, e.g. an exponential distribution? (Or for whichever distributions it is relevant.)
Sometimes? As far as I understand, distributions.jl simply refuses to deal in Float32. I would accept “it’s good enough, it’s impossible to distinguish in Float64 by sampling”. Using 40 bits of randomness instead of 52 would also be OK by the same standard.
Stdlib already provides `Random.randexp`, so `log(rand())` is already an RTFM-failure anyways.
But sure, somebody ought to audit that package. It’s probably a good canary in the “experts are already aware of this” coalmine.
---
### Author:  [@mikmoore](https://discourse.julialang.org/u/mikmoore)
#### Post date: [October 19, 2023, 3:38pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/11 "2023-10-19T15:38:25Z")
One question is how much entropy we intend to generate these numbers from? Do we generate a `rand(Float32)` from a `rand(UInt32)`, a `rand(UInt64)`, a system `rand(UInt)` (system-specific seems awful), or something bigger?
Let’s consider the case of `Float32`. I might be off by a bit in a spot or three here but forgive me that. We have 23 bits of mantissa. If we attempt to generate this from a `UInt32`, we have 9 “spare” bits we can use to extend the dynamic range beyond just `rand12()-1` (or `rand(0:2^23-1) .* 2^-23`). “Unfortunately,” those 9 bits are uniformly distributed but our exponent bits are not. Each smaller octave is half as likely to be utilized as the one above it. So really, those 9 bits means we can extend our range an extra 9 octaves (certainly not 2^9=512 octaves).
An easy way to turn those top 9 bits into an appropriately distributed exponent would be to use `leading_zeros`/`leading_ones`, since they produces 0 with 50% probability, 1 with 25%, etc. We’ll need to clamp the result of `leading_X` so that it does not read into the significand bits.
That said, isn’t it possible to just use a native machine conversion of a `UIntXX` to a `FloatXX` (except for the 16 bit case, since it can overflow), then scale by an appropriate power of two to get a value in [0,1)? The hardware should truncate the bits properly and maintain the distribution, no?
---
### Author:  [@mbauman](https://discourse.julialang.org/u/mbauman)
#### Post date: [October 19, 2023, 3:41pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/12 "2023-10-19T15:41:46Z")
I mean, it’s all a matter of tradeoffs. As I see it, there are really three key questions that’d need to be considered:
- Is it breaking? It changes the interval from closed-open to open-open, and it’s a change in the documented behaviors. We have said, though, that changing RNG results is ok. It’s worth reviewing the time we looked into unbiasing `rand()` by making _just_ the change from half-open to open: [make rand unbiased for floats #33222](https://github.com/JuliaLang/julia/issues/33222). My read there is that _maybe_ but the weighting of tradeoffs in that change (which is different than this) weren’t worth it.
- Complexity and correctness. It’s quite clear what the distribution of `(0:2^53-1) ./ 2^53` looks like, flaws and all — and there’s additionally very well-established literature and norms in other languages. Apparently zig didn’t _quite_ get its newfangled implementation right: [Oh, cute. It looks like they ever-so-slightly overweight the probability of valu... | Hacker News](https://news.ycombinator.com/item?id=37934946).
- Performance.
---
### Author:  [@adamslc](https://discourse.julialang.org/u/adamslc)
#### Post date: [October 19, 2023, 4:45pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/13 "2023-10-19T16:45:22Z")
Even if there is no change to code, I do think that an edited version of this post deserves a spot in the manual. I had no idea about these problems, and this post was about the perfect level of complexity to educate me.
---
### Author:  [@danielwe](https://discourse.julialang.org/u/danielwe)
#### Post date: [October 19, 2023, 5:45pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/14 "2023-10-19T17:45:57Z")
> [@Oscar\_Smith](#):
>
> the multiply approach gets 1 extra bit of randomness
For this reason, I made a suggestion a while back to make this conversion the fallback option for other RNGs as well: [Fallback `rand(Float32/64)` discards more entropy than necessary · Issue #44887 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/issues/44887). But a replacement that increases entropy even more and improves rare event statistics in the small floats would of course be even better.
---
### Author:  [@Sukera](https://discourse.julialang.org/u/Sukera)
#### Post date: [October 19, 2023, 6:09pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/15 "2023-10-19T18:09:25Z")
> [@mbauman](#):
>
> Is it breaking? It changes the interval from closed-open to open-open, and it’s a change in the documented behaviors. We have said, though, that changing RNG results is ok.
We guarantee the distribution, so no, we cannot remove `0`.
See also [here](https://docs.julialang.org/en/v1/stdlib/Random/#Reproducibility):
> The statistical _distribution_ from which random samples are drawn _is_ guaranteed to be the same across any minor Julia releases.
---
### Author:  [@mkitti](https://discourse.julialang.org/u/mkitti)
#### Post date: [October 19, 2023, 6:19pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/16 "2023-10-19T18:19:10Z")
Could we introduce this as a new random number generator?
---
### Author:  [@danielwe](https://discourse.julialang.org/u/danielwe)
#### Post date: [October 19, 2023, 6:20pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/17 "2023-10-19T18:20:59Z")
Thanks for bringing this up. I discovered this rabbit hole a while ago while debugging an issue, as documented here: [Overflow in `randn` using CUDA.jl's native RNG · Issue #1464 · JuliaGPU/CUDA.jl · GitHub](https://github.com/JuliaGPU/CUDA.jl/issues/1464#issuecomment-1090893798), but didn’t take it any further and I am glad that someone more knowledgeable is on it.
Just as a point of comparison, it’s clear based on my experiments that Nvidia’s cuRAND API is doing something smarter than the simplest conversions, at least for Float32, but it’s not going quite as deep as the algorithm proposed on HN. Specifically, based on the smallest value I could observe being 2^{-33} and its recurrence rate about the same, my guess is that it takes 32 random bits `r::UInt32` and rounds (2r + 1) / 2^{33} to the nearest representable float without any intermediate rounding.
Note also that cuRAND is explicitly documented to return values in (0, 1]:
> The curandGenerateUniform() function is used to generate uniformly distributed floating point values between 0.0 and 1.0, where 0.0 is excluded and 1.0 is included.
These discrepancies are perhaps relevant to the discussion of how much people rely on the low-level specifics of `rand` and how breaking a change would be. Clearly, commonly used APIs for different hardware choose different conventions.
(Note if you’re doing your own experiments: CUDA.jl does not use cuRAND by default, but makes it available through `CUDA.CURAND.default_rng()`.)
---
### Author:  [@danielwe](https://discourse.julialang.org/u/danielwe)
#### Post date: [October 19, 2023, 6:29pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/18 "2023-10-19T18:29:26Z")
> [@Sukera](#):
>
> We guarantee the distribution, so no, we cannot remove `0`.
>
> See also [here](https://docs.julialang.org/en/v1/stdlib/Random/#Reproducibility):
>
> > The statistical _distribution_ from which random samples are drawn _is_ guaranteed to be the same across any minor Julia releases.
If you interpret “statistical distribution” here as specifying the exact Float32/Float64 values on which the distribution has support, this was already broken when the Xoshiro RNG was introduced in 1.7 with the slightly improved bits-\>float conversion whose support is twice as dense:
> [@rfourquet](#):
>
> Actually, specifically for `Xoshiro`, we generate floats slightly differently
---
### Author:  [@Sukera](https://discourse.julialang.org/u/Sukera)
#### Post date: [October 19, 2023, 6:39pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/19 "2023-10-19T18:39:59Z")
> [@danielwe](#):
>
> If you interpret “statistical distribution” here as specifying the exact Float32/Float64 values on which the distribution has support, this was already broken when the Xoshiro RNG was introduced in 1.7
I think it’s meant as the `[0, 1)` interval, not the exact output set. That is, it’s meant as an interval on the real number line, not on floating points.
---
### Author:  [@danielwe](https://discourse.julialang.org/u/danielwe)
#### Post date: [October 19, 2023, 6:51pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/20 "2023-10-19T18:51:19Z")
That makes sense. But as explained in the OP, if you exploit the full density of small floats (even excluding subnormals) you would never actually observe zero even when sampling over [0, 1)—the subset of this interval that would round to zero rather than to a finite float would be so vanishingly small. But you could of course maintain a practically dead code path returning zero to formally preserve backward compatibility.
[Next page](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184.md?page=2)