# 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:** 2

<div class="post-metadata">

**Author:** ![Sukera](https://avatars.discourse-cdn.com/v4/letter/s/ce7236/32.png) [@Sukera](https://discourse.julialang.org/u/Sukera)\
**Post date:** [October 19, 2023, 8:01pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/21 "2023-10-19T20:01:33Z")

</div>

Sure, but at the same time you can say the same thing of any floating point number that we currently can generate. Or of any `Int64`, for that matter; the vast majority of them will be large, simply because the upper bits have more “weight” to them. That doesn’t mean that it’s ok to just remove any one number completely though; you’re changing the distribution. Consider changing the example from the OP, only testing for `0.5f0` instead of `0.0`:

```julia
julia> function func()
              N = 0
              while true
              N += 1
              rand(Float32) == 0.5f0 && return N
              end
              end
func (generic function with 1 method)

julia> func()
1532588

julia> func()
4351377

julia> func()
39273158

```

Curious, isn’t it? The documented promise we make is that these numbers are all (roughly) in the same ballpark, no matter which `x` you choose to compare to. The example works for any `Float32` outside of subnormals (because we don’t generate those). Whether that distribution is uniform or not cannot be checked by this kind of spot-checking though - you actually have to go in and create a histogram of the numbers; shouldn’t be too difficult, since that’s only about 17GB of data (32 bits per Float32 means ~4.3 billion numbers, if you store the counts in a `UInt32` you end up with ~17GB of data). The claim in the OP that this is a “correctness footgun” is so far (to my eyes at least) unfounded.

So to me at least, the question posed in the OP is less about “should we generate 0.0f0” or “is the default `rand` uniform” and more about “should we generate subnormals” - and for that I say “not in the default, but we ought to provide a generator specifically for those” (maybe in a package?), because once you do, you very easily get a dramatic performance hit due to bad hardware support for them. This can also be a security problem [[1](https://cseweb.ucsd.edu/~dkohlbre/papers/subnormal.pdf)].

Additionally, I don’t think we ever guaranteed to generate every possible `Float32` or `Float64` in the real interval `[0, 1)`, and I really wouldn’t expect us to by default, primarily for these performance reasons down the line (and not because generating the numbers itself is a tad bit slower).

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [October 19, 2023, 8:18pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/22 "2023-10-19T20:18:14Z")

</div>

> [@Sukera](#):
>
> The documented promise we make is that these numbers are all (roughly) in the same ballpark, no matter which `x` you choose to compare to.

If I’m understanding you here that’s not the case, and no, this isn’t about subnormals at all as neither scheme generates them. Your `func` will never terminate, for example, if you look for `prevfloat(0.5f0)`.

It would exist in the output (but be **half as likely to occur** ) in the scheme proposed by @foobar_lv2 and the above blog post. But this makes some sense if you think about how floats represent rounded “bins” of real numbers, and the bins in [0.25, 0.5) are half as big.

---

<div class="post-metadata">

**Author:** ![Sukera](https://avatars.discourse-cdn.com/v4/letter/s/ce7236/32.png) [@Sukera](https://discourse.julialang.org/u/Sukera)\
**Post date:** [October 19, 2023, 8:30pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/23 "2023-10-19T20:30:28Z")

</div>

> [@mbauman](#):
>
> If I’m understanding you here that’s not the case, and no, this isn’t about subnormals at all as neither scheme generates them. Your `func` will never terminate, for example, if you look for `prevfloat(0.5f0)`.

That’s a good point, and would have been better as an example in the OP than `0.0f0`. Regardless, my point stands - we haven’t guaranteed to generate every possible floating point number (neither normal nor subnormal, I don’t think) in the real interval `[0, 1)`.

> [@mbauman](#):
>
> It would exist in the output (but be **half as likely to occur** ) in the scheme proposed by @foobar_lv2 and the above blog post. But this makes some sense if you think about how floats represent rounded “bins” of real numbers, and the bins in [0.25, 0.5)[0.25,0.5)[0.25, 0.5) are half as big.

That’s true - but changing that to be bin-aware and actually generating every (normal) number possible in the bin definitely changes the distribution (as I understand it), and that’s not something we can do. The guarantee is on the distribution of the real line, not the floating point output, because that’s the distribution actually relevant when thinking about statistical guarantees of the desired randomness.

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [October 19, 2023, 9:30pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/24 "2023-10-19T21:30:39Z")

</div>

I don’t understand your point, then. Those last two sentences there really seem to be at odds with each other, no? If the guarantee is on the distribution of the real line and not the floating point output, then the likelihood of an exact value being generated isn’t guaranteed. And, if you take the perspective of all the real numbers from [0, 1), the chances of a random one being closer to `0f0` than any other nonzero value are on the order of 10^{-45} — astronomically unlikely. So are we really changing the distribution? Or are we making it a more faithful approximation the uniform real numbers?

Admittedly I agree, the [blog post](https://www.corsix.org/content/higher-quality-random-floats) and exposition here weren’t the easiest to follow, and I definitely wouldn’t call the current generator the result of catastrophic cancellation. But the core idea is quite straightforward:

- The status quo is that we (effectively) sample from `0:eps(T)/2:1-eps(T)/2` with equal probability.
- This change would sample from `0.5:eps(0.5):1-eps(0.5)` with probability \omega, `0.25:eps(0.25):0.5-eps(0.25)` with probability \frac{\omega}{2}, `0.125:eps(0.125):0.25-eps(0.125)` with probability \frac{\omega}{4}, and so on… with 11 \* I think such regions for Float64 and 8\* for Float32. This means that **all values** greater than `1.3234889800848443e-23` and `1.1641532f-10` are possible outputs.

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [October 19, 2023, 9:40pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/25 "2023-10-19T21:40:43Z")

</div>

On a slight tangent, `Uniform(a,b)` seems to need a documentation fix, as it says:

```julia
Uniform(a,b)

  The continuous uniform distribution over an interval [a, b] has probability density function

  f(x; a, b) = \frac{1}{b - a}, \quad a \le x \le b

  Uniform() # Uniform distribution over [0, 1]
  Uniform(a, b) # Uniform distribution over [a, b]

```

yet, it uses `a + rand()*(b-a)` which is kind of `[a, b)`.

---

<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, 10:00pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/26 "2023-10-19T22:00:13Z")

</div>

> [@mbauman](#):
>
> and 8\* for Float32

125 such regions if we stick to non-subnormal Float32.

> [@mbauman](#):
>
> Admittedly I agree, the [blog post](https://www.corsix.org/content/higher-quality-random-floats) and exposition here weren’t the easiest to follow,

I’d be super happy if you could find better words to explain the issue, or maybe describe what makes it hard to follow?

What is the conceptual difficulty that makes this hard to grok?

Are people just unaware how floating point numbers work, conceptually? That’s why I added the reminder on floating point numbers, but I admit that my post will be much too concise for people who need an intro-to-float instead of a reminder-of-float.

---

<div class="post-metadata">

**Author:** ![adienes](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/adienes/32/37459_2.png) [@adienes](https://discourse.julialang.org/u/adienes)\
**Post date:** [October 19, 2023, 10:08pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/27 "2023-10-19T22:08:14Z")

</div>

— not saying these are the options, just setting up a strawman —

if a change to `rand()` were proposed that made the distribution more uniform, or with more bits of randomness, or more performant

I don’t think changing the _precise_ categorical distribution over the set of floats is the kind of “breaking” change that will get (m)any complaints, besides those cases where known bias corrections are now over-correcting

---

<div class="post-metadata">

**Author:** ![Sukera](https://avatars.discourse-cdn.com/v4/letter/s/ce7236/32.png) [@Sukera](https://discourse.julialang.org/u/Sukera)\
**Post date:** [October 19, 2023, 10:13pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/28 "2023-10-19T22:13:43Z")

</div>

> [@mbauman](#):
>
> I don’t understand your point, then. Those last two sentences there really seem to be at odds with each other, no?

Hm, maybe this can illustrate my point; the bin between 1-2 we currently use is a whole bin, there are no smaller subdivisions. That is, the distance between each successive number in that bin is equal. We effectively then “shift” that bin towards zero on the real number line by subtracting one, keeping the distance between numbers the same. With the approach from the article, that doesn’t happen anymore; due to truly picking from the available _floating point_ numbers, the numbers end up denser the closer you are to zero. There simply are more numbers available, and that’s not necessarily a good thing.

The numbers are still all equally likely to be picked (weighted by the size of the rounding bins they represent), but if you actually plot the numbers you get from that random generator, you’ll notice a cluster at the bottom (irrespective of higher bins being more likely and those going to be chosen more often! In a histogram, all numbers from the possible set of _floating point_ numbers are perfectly equally likely). That is the shift in distribution I’m talking about.

I may very well have something wrong in my logic, and I admittedly haven’t plotted the results (which would be easiest way to convince people either way), but this is what I’m talking about.

* * *

Another way to put it - keeping the current local density around any given number we generate the same is IMO more important than being able to generate all possible _floating point_ values.

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [October 19, 2023, 10:18pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/29 "2023-10-19T22:18:15Z")

</div>

> [@Sukera](#):
>
> but changing that to be bin-aware and actually generating every (normal) number possible in the bin definitely changes the distribution (as I understand it)

This seems like a misunderstanding. Nobody is suggesting that each floating point value in [0, 1) should have equal weight. Rather, each value in sample space should be weighted by the width of the bin it represents on the real line. This width is determined by the distance to the nearest neighbors in sample space.

> [@Sukera](#):
>
> we haven’t guaranteed to generate every possible floating point number (neither normal nor subnormal, I don’t think) in the real interval `[0, 1)`.

Nope, but the OP explains why it would be good to do so, or at least to improve on the current situation. It probably doesn’t matter in cases where U([0, 1)) is the actual distribution you’re interested in (e.g., for the `if rand() < p` trick), but it makes a huge difference for rare event statistics if you’re actually trying to sample a heavy-tailed distribution by feeding your uniform rand through the inverse CDF.

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [October 19, 2023, 10:23pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/30 "2023-10-19T22:23:56Z")

</div>

> [@Sukera](#):
>
> you’ll notice a cluster at the bottom

Actually, if you sample enough values and study them closely enough, you’ll notice that they are more clustered in the top end, not the bottom, because the sample space is sparser up there. But any histogram with bin size larger than the largest rounding bin will be flat (given enough samples).

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [October 19, 2023, 10:24pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/31 "2023-10-19T22:24:23Z")

</div>

> [@foobar\_lv2](#):
>
> 125 such regions if we stick to non-subnormal Float32.

Ah, yes, of course it couldn’t be 8 regions (as that’d have a minimum value of 2^{-8}). But it’s also not 125 (at least not in your implementation that I played with yesterday). This is representative of my own confusion here — things are a bit of a mishmash of bit twiddling with multiple implementations, floating point shenanigans, passionate pushes about good and bad, and theoretical probabilities.

On the bit twiddling front: with 23 bits for the significand, there are only 8 or 9 bits left to use for the exponent, no? It appears you’re using all 32 bits in the lzcnt (the minimum value is 2^{-33})… which would impart quite the bias… so I’m obviously missing something here.

---

<div class="post-metadata">

**Author:** ![Sukera](https://avatars.discourse-cdn.com/v4/letter/s/ce7236/32.png) [@Sukera](https://discourse.julialang.org/u/Sukera)\
**Post date:** [October 19, 2023, 10:26pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/32 "2023-10-19T22:26:16Z")

</div>

> [@danielwe](#):
>
> This seems like a misunderstanding. Nobody is suggesting that each floating point value in [0, 1) should have equal weight. Rather, each value in sample space should be weighted by the width of the bin it represents on the real line. This width is determined by the distance to the nearest neighbors in sample space.

Yes, I’m aware. I’m not talking about that though.

> [@danielwe](#):
>
> Actually, if you sample enough values and study them closely enough, you’ll notice that they are more clustered in the top end, not the bottom, because the sample space is sparser up there. But any histogram with bin size larger than the largest rounding bin will be flat (given enough samples).

I’m not talking about a histogram. Yes, if you plot a histogram there will be a cluster at the top (which is also bad). I’m talking about “did a number occur there at all” which has a cluster close to zero. That is not a good clustering, because it indicates that the numbers closer to 1 we can actually represent are less diverse than the numbers close to 0 - i.e., they don’t represent a uniform sampling, but a biased one.

You can also see exactly that in the histogram you describe.

* * *

Ultimately, the question here boils down to “sample from the _floating point_ values representable in `[0, 1)`” or “sample from the real number line with uniform distances (== bin sizes) in `[0, 1)`”. My point is that we currently assume/guarantee the latter, while the article and OP advocate for the former.

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [October 19, 2023, 10:38pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/33 "2023-10-19T22:38:40Z")

</div>

> [@Sukera](#):
>
> I think it’s meant as the `[0, 1)` interval,

This has the same measure as (0,1)

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [October 19, 2023, 10:50pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/34 "2023-10-19T22:50:26Z")

</div>

Floats are logarithmically rather than linearly spaced because this is the useful and appropriate representation for arithmetic across orders of magnitude. It’s the reason you can do `1/x` with impunity. This is still true if `x` originated from an RNG. In general, increasing the entropy of random numbers can only be a good thing, apart from performance regressions and [xkcd: Workflow](http://xkcd.com/1172) issues. Of course, any proposed conversion algorithm should be carefully analyzed for statistical bias.

* * *

Btw., aren’t these interpretations contradictory?

> [@Sukera](#):
>
> I think it’s meant as the `[0, 1)` interval, not the exact output set.

> [@Sukera](#):
>
> My point is that we currently assume/guarantee […] “sample from the real number line with uniform distances (== bin sizes) in `[0, 1)`”

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [October 19, 2023, 10:51pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/35 "2023-10-19T22:51:57Z")

</div>

Other thoughts. If X is uniform on (0,1) and Y is Normal(0, 1e-6) then X+Y restricted to (0,1) (ie. Rejection sampling) is also uniform on (0,1) but it’s more densely distributed on the floats.

---

<div class="post-metadata">

**Author:** ![Sukera](https://avatars.discourse-cdn.com/v4/letter/s/ce7236/32.png) [@Sukera](https://discourse.julialang.org/u/Sukera)\
**Post date:** [October 19, 2023, 11:03pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/36 "2023-10-19T23:03:35Z")

</div>

> [@danielwe](#):
>
> Floats are logarithmically rather than linearly spaced because this is the useful and appropriate representation for arithmetic across orders of magnitude.

Sure, but again, that then changes the question of what is sampled from. Are we sampling from the representable floating point values in `[0, 1)` (i.e., picking one of a logarithmically sized bin), or are we sampling (mathematically) from the real number line and binning to uniform bin sizes afterwards? I think the “uniform distribution” in the docs is talking squarely about the latter, and that should be represented in the binning of the actual values we generate.

If you take the sampled-from-representable-floating-point-values distribution and re-bin it to have uniform bin sizes you of course get back a nice flat histogram - but then you could have just generated that in the first place by using the existing method. Not doing that on the other hand skews your histogram due to the uneven binning, and some actually representable values are then more likely to occur than others 🤷

> [@danielwe](#):
>
> Of course, any proposed conversion algorithm should be carefully analyzed for statistical bias.

Yes, and I’m worried about bias being introduced due to the difference in bin sizes used across the sampling domain.

> [@danielwe](#):
>
> Btw., aren’t these interpretations contradictory?

Good catch - the former certainly doesn’t imply the latter (though I had the latter in mind while writing the former - sorry about that!).

Again though, it comes down to where & how to do the binning. I think sticking to uniform bin sizes is a better default because it keeps the output likelihood of any single floating point value we are able to generate constant. It’s just one less potential source of bias.

---

<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 20, 2023, 10:25am UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/37 "2023-10-20T10:25:31Z")

</div>

> [@danielwe](#):
>
> (e.g., for the `if rand() < p` trick)

That’s actually a good one. Consider

```julia
import Random
struct Bernoulli{T<:AbstractFloat}
p::T
end
sample(b::Bernoulli) = sample(Random.default_rng() ,b)
sample(rng, b::Bernoulli{T}) where T = rand(rng, T) < b.p

myRareCoin = Bernoulli(1f-8)

```

Does this code scream NOPE NOPE NOPE to you @Sukera ?

```julia
julia> function testRand(b::Bernoulli, N)
       hits = 0
       for i in 1:N
       sample(b) && ( hits += 1 )
       end
       @show b, hits, N, hits/N
       end

julia> @time testRand(Bernoulli(1f-8), 1<<32)
(b, hits, N, hits / N) = (Bernoulli{Float32}(1.0f-8), 265, 4294967296, 6.170012056827545e-8)
  5.170603 seconds (58 allocations: 3.164 KiB)
(Bernoulli{Float32}(1.0f-8), 265, 4294967296, 6.170012056827545e-8)

julia> @time testRand(Bernoulli(1e-8), 1<<32)
(b, hits, N, hits / N) = (Bernoulli{Float64}(1.0e-8), 44, 4294967296, 1.0244548320770264e-8)
  5.095529 seconds (58 allocations: 3.453 KiB)
(Bernoulli{Float64}(1.0e-8), 44, 4294967296, 1.0244548320770264e-8)

```

My rare coin is overweight by a factor of 6.

---

<div class="post-metadata">

**Author:** ![Sukera](https://avatars.discourse-cdn.com/v4/letter/s/ce7236/32.png) [@Sukera](https://discourse.julialang.org/u/Sukera)\
**Post date:** [October 20, 2023, 2:06pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/38 "2023-10-20T14:06:16Z")

</div>

> [@foobar\_lv2](#):
>
> Does this code scream NOPE NOPE NOPE to you @Sukera ?

Considering the current generation scheme more or less has

```julia
julia> eps(nextfloat(1.0f0))
1.1920929f-7

julia> eps(prevfloat(2.0f0))
1.1920929f-7

```

as its accuracy, yes. I agree though that requiring to know this is neither ergonomic nor good for a user 🤔 On the other hand, one could also argue that this is inevitable depending on which value you pick; don’t you get the same “overweighing” on the other end of the interval, no matter whether this is fixed for small numbers? After all, the “proper” bin has a higher precision than the interval we currently generate too:

```julia
julia> eps(prevfloat(1f0))
5.9604645f-8

```

Either way, if this is changed, it should definitely be clarified that `rand()` samples from the representable floating point values in the interval `[0, 1)`, and not from the real number line.

---

<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 20, 2023, 2:32pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/39 "2023-10-20T14:32:22Z")

</div>

> [@Sukera](#):
>
> On the other hand, one could also argue that this is inevitable depending on which value you pick; don’t you get the same “overweighing” on the other end of the interval, no matter whether this is fixed for small numbers? After all, the “proper” bin has a higher precision than the interval we currently generate too:

This is not inevitable.

A user who writes `Bernoulli(1.0f0 - 1.0f-8)` or more generally `Bernoulli(p)` for p close to 1 is very hard to help – on some level this is a question of basic computer literacy and the user “deserves” what’s coming to them.

My example (`p` close to 0) bites a user who does everything correctly. The expected deviation for `Bernoulli(1.0f-8)` is on order `eps(1.0f-8) == 8.881784f-16`.

In other words, the implementation `sample(rng, b::Bernoulli{T}) where T = rand(rng, T) < b.p` is the bug. This implementation looks intuitively right, also when considering all the documentation of `rand`. Nevertheless, this implementation is incorrect, because either the implementation or documentation of `rand` is bad.

Btw @adienes , Rust has a similar issue: There is no mismatch between docs and impl in rust-rand; but the docs are not energetic enough about warning users about the lacking dynamic range. I mean, not even rust-rand itself understands the (documented!) problem with its uniform rand numbers, to the point that pareto misuses it and creates completely wrong (fat) tails.

---

<div class="post-metadata">

**Author:** ![mikmoore](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikmoore/32/31109_2.png) [@mikmoore](https://discourse.julialang.org/u/mikmoore)\
**Post date:** [October 23, 2023, 4:49pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/40 "2023-10-23T16:49:33Z")

</div>

Suppose we sample a float v \in [0,1). Define the smallest nonzero value we may observe as \eta. Afterward, we dutifully verify that v \notin \{0\} so that v \in (0,1). We might now reasonably expect that 1-v \in (0,1). But if 1-\eta rounds to 1 (e.g., because we apply “extra” precision to small v), we instead get 1-v \in (0,1]. Obviously, this is now the domain of floating point pitfalls more than random number generation, but it’s a sharp edge that a casual user might reasonably trip on. I imagine there are a number of related inconsistencies, possibly more insidious than this.

It’s usually not difficult to construct failure modes with floating point math, and a careless implementation will often run afoul at least one. My example doesn’t seem _intrinsically_ more or less an issue than the “generate small-probability Bernoulli random variables” issue raised above. I don’t see how the two can be simultaneously resolved by any adjustment to the set of values returned by `rand`. Neither resolution avoids pitfalls for all possible uses.

* * *

This reminds me of a snippet of code I had in a recent project. I needed to compute \log(1-\exp(x)) for x\le0. My implementation was `x > -1 ? log(-expm1(x)) : log1p(-exp(x))`. Note that the first calculation is unsuitable when x \ll -1 and the second is unsuitable when x \approx 0. No single code path was suitable for the entire domain (at least not with standalone `log`/`exp`-like functions), so extra care was necessary.

If you truly want accurate `Bernoulli(1f-40)` random variables, applying a few bits of extra precision won’t fix a `rand() < p` implementation. The sampler itself must find a way to produce a reasonable value (or error) for such a challenging request. This probably means switching to a different algorithm in certain regimes.

By the reasoning I’ve discussed here, and status-quo bias, I’m marginally more sympathetic to the current “fixed-point” output precision. Though I’m willing shift my stance given further examples.

[Previous page](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184.md?page=1)

[Next page](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184.md?page=3)