# 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:** 1\
**Showing post:** 1

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@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")

</div>

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!

---

_[View the full topic](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184)._