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

<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 29, 2023, 2:28am UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/222 "2023-10-29T02:28:07Z")

</div>

> [@dlakelan](#):
>
> so I still like the 0 fixup.

Why not just do `b = (i >> 32) + 1`. Adds a 1/2^64 glitch around 1/2^32 but saves a conditional (maybe not so crucial for speed - needs measurement).

---

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [October 29, 2023, 4:54am UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/223 "2023-10-29T04:54:58Z")

</div>

> [@foobar\_lv2](#):
>
> So I am arguing that, if we take 64 bits of random due to performance tradeoffs, the smallest achievable number should be of order `ldexp(1.0, -64)`, in order to avoid the extreme overweighting of small numbers.

This is a good line of argument that I want to continue. Not in a position to do so right now but this feels like the right direction to discuss.

---

<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 29, 2023, 5:33am UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/224 "2023-10-29T05:33:19Z")

</div>

@StefanKarpinski I agree, and it’s more or less the argument I thought that was always being made, so I like that we’re converging on some common understanding.

Here’s some data on timing for my slightly modified code above:

```julia
julia> function rand32()
                i = rand(UInt64)
                a = i & 0xffffffff
                b = i >> 32
                if a == 0x0
                  if b == 0 ; b= 1 ; end
                  c = ldexp(Float32(b),-64)
                 else 
                   c = ldexp(Float32(a),-32)
                end
                return c
             end
rand32 (generic function with 1 method)

julia> @btime(rand32())
  4.308 ns (0 allocations: 0 bytes)
0.4268289f0

julia> @btime(rand(UInt64))
  2.866 ns (0 allocations: 0 bytes)
0xbcb048922324e5f6

julia> 4.308/2.866
1.5031402651779482

```

So it’s about 1.5 times the cost compared to just generating the 64 bits, and improves the left tail a LOT, the smallest nonzero value it can generate is ldexp(1.0f0,-64) compared to I think ldexp(1.0,-24) for the current code?

The smallest Float32 on the other hand is about 2.58f-26 as big as ldexp(1.0f0,-64)

So we have the option also to extend the dynamic range even further, we could do something like ldexp(Float32(rand(UInt32)),-64-32) as a 3rd step of extending the dynamic range even farther, it would only get hit once every 2^64 calls.

---

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

</div>

How about:

```julia
function rand32_while()
    exponent = -32
    while ((a = rand(UInt32)) == 0)
        exponent -= 32
    end
    return ldexp(Float32(a), exponent)
end

```

(chance of `0.0` is around `2^(-128)`)

or same semantics but simpler to prove finite run-time:

```julia
function rand32_for()
    local exponent, a
    for outer exponent in -32:-32:-196
        (a = rand(UInt32)) == 0 || break
    end
    return ldexp(Float32(a), exponent)
end

```

---

<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 29, 2023, 3:02pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/226 "2023-10-29T15:02:27Z")

</div>

> [@Dan](#):
>
> `rand(UInt32)`

Just from a speed perspective, my understanding was the RNG generates 64 bits at a time, so I don’t know if we are throwing those away and calling the generator more times than needed? Maybe it stores the unused bits and returns them quickly at the second call?

Note, maybe it doesn’t matter, we only make a second call once every 4 billion (2^32) calls, and a third call once every 2^64

---

<div class="post-metadata">

**Author:** ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)\
**Post date:** [October 29, 2023, 3:13pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/227 "2023-10-29T15:13:14Z")

</div>

> [@dlakelan](#):
>
> we only make a second call once every 4 billion (2^32) calls

Presumably a real implementation would eliminate the branching/loop anyway by getting all the random bits upfront?

---

<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 29, 2023, 3:17pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/228 "2023-10-29T15:17:49Z")

</div>

> [@nsajko](#):
>
> getting all the random bits upfront?

Which is better, avoiding a branch but calling the RNG twice or 3 times as much as needed every time, or avoiding the second RNG almost always? I guess time to benchmark. I’m guessing the loop is probably faster since it almost always goes once and the RNG won’t SIMD per previous post.

---

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

</div>

```julia
julia> @btime(rand32_for())
  4.067 ns (0 allocations: 0 bytes)
0.15846866f0

julia> @btime(rand(UInt64))
  2.875 ns (0 allocations: 0 bytes)
0x7ae0be730a4bf57b

julia> @btime(rand(UInt32))
  2.865 ns (0 allocations: 0 bytes)
0x43bc04e7

```

So, I love it. I also like the idea of doing the same thing for Float64 since it shouldn’t be the case that Float64 is more problematic than Float32

```julia
julia> function rand64_for()
           local exponent, a
           for outer exponent in -64:-64:-256
               (a = rand(UInt64)) == 0 || break
           end
           return ldexp(Float64(a), exponent)
       end
rand64_for (generic function with 1 method)

julia> @btime(rand64_for())
  4.548 ns (0 allocations: 0 bytes)
0.26893953364756196

julia> @btime(rand(Float64))
  2.845 ns (0 allocations: 0 bytes)
0.8204255370807716

julia> 4.548/2.845
1.5985940246045693

```

---

<div class="post-metadata">

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

</div>

> [@dlakelan](#):
>
> `@btime(rand32_for())`

This tries to measure the latency of a single call, but I guess it could be more interesting to benchmark applying `rand` to an array.

---

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

</div>

> [@nsajko](#):
>
> This tries to measure the latency of a single call, but I guess it could be more interesting to benchmark applying `rand` to an array.

There are really 3 distinct cases we care about:

1. Calling `rand` often, but in between there is something that the compiler cannot optimize through. In this case, the internal state of the RNG is read, modified, and written back for each call. The read/write probably doesn’t hit L1 cache and instead goes via store buffer.
2. Calling `rand` in a tight loop that the compiler understands. The internal RNG state is read only once and written only once, otherwise it is kept in register. Register is faster than store buffer.
3. Calling `rand!` on an array. This is a completely separate codepath and algorithm that uses SIMD.

---

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

</div>

So, to summarize, we have a simple candidate algorithm for rand() for Float64 and Float32 which uses a finite loop, is guaranteed to terminate, and generates 0.0 at less than 2^-128 = 2.9e-39 for either type, which means a million core billion hz supercomputer generates it less than once every 10^16 years in the absence of software bugs. It also has the CDF property that cdf(x) = x for values between 0.0 and 1.0 to within the capacity of the RNG to be fully uniform, and it costs about 1.5 times as much as generating a single UInt64.

I consider this a job well done. What’s the consensus?

---

<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 29, 2023, 8:24pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/233 "2023-10-29T20:24:29Z")

</div>

> [@dlakelan](#):
>
> I consider this a job well done. What’s the consensus?

Pretty good.

Two points:

1. There is a related conundrum I’ve encountered:

```julia
julia> @benchmark Float32(rand(UInt64))
BenchmarkTools.Trial: 10000 samples with 1000 evaluations.
 Range (min … max): 6.384 ns … 30.382 ns ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 6.842 ns ┊ GC (median): 0.00%
 Time (mean ± σ): 6.878 ns ± 0.384 ns ┊ GC (mean ± σ): 0.00% ± 0.00%

           ▂▅▆██▇▆▃                                           
  ▂▂▂▂▃▃▄▆██████████▇▄▄▃▂▂▂▂▂▂▁▂▁▂▂▂▂▁▁▂▂▂▁▁▂▂▂▂▁▂▂▂▂▂▂▂▂▂▂▂ ▃
  6.38 ns Histogram: frequency by time 8.39 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

```

```julia
julia> @benchmark Float32(rand(UInt32))
BenchmarkTools.Trial: 10000 samples with 1000 evaluations.
 Range (min … max): 3.075 ns … 150.727 ns ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 3.159 ns ┊ GC (median): 0.00%
 Time (mean ± σ): 3.190 ns ± 1.488 ns ┊ GC (mean ± σ): 0.00% ± 0.00%

                           ▁ ▁▆▂▆▃▄█▂▇▂▃▇▃▃▇▁▄                 
  ▂▁▁▁▁▂▂▁▁▂▂▂▂▂▂▂▂▂▂▃▃▃▄▆▅█▇█████████████████▇▅▅▄▄▄▃▄▃▃▂▂▂▂▂ ▄
  3.08 ns Histogram: frequency by time 3.22 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

```

The two are identical in terms of `rand` calculation, so why does conversion cost more for UInt64 ? Looking at the generated code, there is another branch for the UInt64 case, but I’m not sure what it is for. Anyone?

1. Officially I’m renaming Float16 as BraindeadFloat16. The Float16 cannot hold 0xFFFF which is the the maximum UInt of the same size. IEEE have made the wrong decisions and BFloat16 should replace it all. My theory for why Float16 is this way, is because it has 10 bits of significant means it reached 1024, which means it can be translated into 3 decimal digits. It looks good for all the decimal system fans. Today humans are in much less regard and I think it wouldn’t have gone the same way.

---

<div class="post-metadata">

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

</div>

> [@Dan](#):
>
> Float16 cannot hold 0xFFFF which is the the maximum UInt of the same size

That’s true for any FP format, except possibly for a trivial FP format which only has a mantissa and is thus actually just an integer…? EDIT: OK, thanks for the clarification.

---

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

</div>

I mean 65535 cannot be represented as a float _with an exponent_! It is Inf !

```julia
julia> Float64(floatmax(Float16))
65504.0

```

```julia
julia> Float16(0xffff)
Inf16

```

or to spell it out: the exponent doesn’t have more than ~`log(#bits)+1` bits.

---

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

</div>

I was also thinking along these lines. However, for reasons discussed earlier, I’d rather throw away some bits than have implicit rounding, and I’d like to avoid branches, so here’s the way I would do it:

```julia
function randf32()
    u = rand(UInt64)
    a = Float32(u >>> 40) * Float32(0x1.0p-24)
    b = Float32((u >>> 16) & 0xffffff) * Float32(0x1.0p-48)
    return ifelse(a === 0f0, b, a)
end

```

This takes the top 24 bits to create a \in [0, 1) and the next 24 bits for b \in [0, 2^{-24}), returning `b` if `a == 0` and `a` otherwise. This could of course be iterated further if desired.

Performance looks alright:

```julia-repl
julia> @btime randf32();
  3.446 ns (0 allocations: 0 bytes)

julia> @btime rand(Float32);
  2.942 ns (0 allocations: 0 bytes)

```

What I don’t like so much about this approach is the step change in density at the chosen cutoff. The sample space remains as coarse as the current implementation down to 2^{-24} and is only refined for numbers smaller than that. I think we should aim to do better but also to avoid rounding and branches.

---

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

</div>

> [@Dan](#):
>
> There is a related conundrum I’ve encountered:

There is no good conversion function uint64-\>float32 in x86, there’s only one for signed integers. It’s why my proposal upthread did `Float32(rand(UInt64)>>1)`.

On the upside, if you do use this, it has a vectorized version on AVX2.

---

<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 30, 2023, 12:41am UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/238 "2023-10-30T00:41:04Z")

</div>

In a scheme like this, you can also stagger the stages such that you don’t exhaust the dynamic range of one stage before switching to the next. Here’s a 3-stage version with no rounding that first samples in [0, 1) with resolution 2^{-24}, but if the outcome is \< 2^{-8} resamples in [0, 2^{-8}) with resolution 2^{-32}, and if _that_ outcome is \< 2^{-16} resamples in [0, 2^{-16}) with resolution 2^{-40}.

```julia
function randf32()
    _24BITS = 0x0000000000ffffff
    _2POW16 = 0x0000000000010000
    u = rand(UInt64)       
    u1 = u >>> 40
    u2 = (u >>> 32) & _24BITS
    u3 = (u >>> 24) & _24BITS
    x = Float32(u3) * Float32(0x1.0p-40)
    x = ifelse(u2 < _2POW16, x, Float32(u2) * Float32(0x1.0p-32))
    x = ifelse(u1 < _2POW16, x, Float32(u1) * Float32(0x1.0p-24))
    return x
end

```

```julia-repl
julia> @btime randf32();
  3.894 ns (0 allocations: 0 bytes)

julia> @btime rand(Float32);
  2.937 ns (0 allocations: 0 bytes)

```

It’s possible to iterate three more times and get down to a resolution of 2^{-64} before exhausting a `UInt64`, but performance suffers since the algorithm takes all branches by design.

In general, the stage overlap can be adjusted to trade resolution near zero for resolution near the transition points.

---

<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 30, 2023, 12:52pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/239 "2023-10-30T12:52:38Z")

</div>

Threadmates, remind me what was wrong with:

```julia
randf32() = abs(ldexp(Float32(rand(Int64)),-63))

```

?

ADDED: Seems to me, the only problem is:

```julia
rand(r::Union{TaskLocalRNG, Xoshiro}, ::SamplerTrivial{CloseOpen01{Float16}}) =
    Float16(Float32(rand(r, UInt16) >>> 5) * Float32(0x1.0p-11))

rand(r::Union{TaskLocalRNG, Xoshiro}, ::SamplerTrivial{CloseOpen01{Float32}}) =
    Float32(rand(r, UInt32) >>> 8) * Float32(0x1.0p-24)

```

in Xoshiro.jl which doesn’t use the full 64-bit to generate smaller floats, which indeed should be fixed (at least in a RNG which generates 64-bit chunks anyway and no tradeoff relevant)

---

<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 30, 2023, 1:09pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/240 "2023-10-30T13:09:41Z")

</div>

> [@Dan](#):
>
> `randf32() = abs(ldexp(Float32(rand(Int64)),-63))`

It rounds-to-even for values of `Int64 > 2^24`, so not only does it return 1.0, even floats in the 2^24 octave over overrepresented by 3x (and each octave thereafter by ((N+2)/N)x.

---

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

</div>

> [@mbauman](#):
>
> It rounds-to-even for values of `Int64 > 2^24`,

But it rounds to nearest first. So which _specific numbers_ are 3x more common? When using 63-bits, this rounding issue is only with numbers ~ 2^(-40) ?

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

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