# 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:** 18
### Author:  [@dlakelan](https://discourse.julialang.org/u/dlakelan)
#### Post date: [November 1, 2023, 7:28pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/343 "2023-11-01T19:28:31Z")
Not really sure either how to do that or that it makes sense, considering this is an empirical cdf based on a sample.
On the other hand, here’s the error in the cdf based on 20\_000 samples for convolution (blue) and rand(Float32) (red)
Another draw based on 40\_000 samples
Just on the basis of a couple of these paths, it doesn’t appear to have any major flaws compared to rand(Float32) across the whole region. It’s CLEARLY much better below 1e-6
This does suggest a kind of statistical test though. Look at the error in the empirical CDF for varying size samples and see if it has appropriate statistics for maximum deviation from 0.
---
### Author:  [@mbauman](https://discourse.julialang.org/u/mbauman)
#### Post date: [November 1, 2023, 7:39pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/344 "2023-11-01T19:39:03Z")
You do realize we’re simultaneously discussing performant algorithms that have exactly 0 error, right?
---
### Author:  [@dlakelan](https://discourse.julialang.org/u/dlakelan)
#### Post date: [November 1, 2023, 7:51pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/345 "2023-11-01T19:51:00Z")
I’m aware @foobar_lv2 proposed some kind of “perfect sampling” algorithm but I was under the impression it was less performant.
Honestly Discourse doesn’t make it super easy to search back through the messages… but I guess it was here and it was faster than I remembered
> [@foobar\_lv2](#):
>
> ```julia
> function idealSample(::Type{T}) where T <: Union{Float64, Float32}
>
> ```
But also, my impression was someone had complained that it wasn’t correct?
Since it’s only about a factor of 2 slower than rand(Float32) maybe I should generate some similar graphs from it…
---
### Author:  [@mbauman](https://discourse.julialang.org/u/mbauman)
#### Post date: [November 1, 2023, 8:06pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/346 "2023-11-01T20:06:24Z")
The challenge isn’t so much scalar `rand` as it is `rand!`. The latter is much harder and where SIMD matters.
We have a pretty-darn-close-to-performant and SIMD-able and proven-to-be-exact equal-bit-width [algorithm](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/261) (it’s about 50% slower for `rand!` on Float32). That puts us at N=32 and N=64 in my table — exactly accurate cdfs for all Float32s ≥ 0.001953125f0 and all Float64s ≥ 0.000244140625.
The new clever bit from just a handful of posts back is if any of the batch of floats generated land outside the densely-supported region (\<1% chance), you just re-run it and scale down to get a finer sampling of the small values from the prior step. And iterate until you’re content about the probabilities (from my N-bits-table).
Still needs to be proven out, but seems quite promising.
---
### Author:  [@mikmoore](https://discourse.julialang.org/u/mikmoore)
#### Post date: [November 1, 2023, 8:27pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/347 "2023-11-01T20:27:13Z")
Below I offer a corrected version of my previous SIMD-able full-precision implementation.
Here is scalar pseudocode
```julia
scale = 2^-32
while scale > 0
u = rand(UInt32)
f = float(u, RoundDown) * scale
u >= 2^24 && return f
scale *= 2^-8
end
f = float(u >>> dropbits) * eps(0f0)
return f
```
and the actual code
```julia
using Random
function trunctofloat(u::U) where U<:Base.BitUnsigned
F = Base.floattype(U)
lomask = typemax(U) >>> (8sizeof(F) - Base.significand_bits(F) - 1) # masking with this number ensures exact representability
hi = F(u & ~lomask)
lo = F(u & lomask)
f0 = hi + lo # combine the parts
roundup = f0 - hi > lo # was the roundoff error of f0 upward?
f = reinterpret(F, reinterpret(U, f0) - roundup) # decrement if rounding was upward
return f
end
block_randf(t::Type,n::Val) = block_randf(Random.default_rng(),t,n)
function block_randf(rng::Random.AbstractRNG, ::Type{T}, N::Val) where {T<:Base.IEEEFloat}
U = Base.uinttype(T) # same-width Unsigned
mineps = eps(zero(T)) # minimum spacing between values that we want to achieve
umax = U(trunctofloat(typemax(U)))
sparebits = trailing_zeros(umax) # how many "extra" bits we have per draw
usebits = 8sizeof(U) - sparebits # how many bits we can possibly consume per draw
redrawbelow = U(exp2(usebits)) # values smaller than this lose precision and must be redrawn
redrawscale = T(exp2(-sparebits)) # if a value fails, scale by this and redraw
scale = T(inv(redrawbelow))
todrop = Int(log2(scale) - log2(mineps)) # how many octaves we need to drop via `redrawscale`
numscales,leftover_bottom = fldmod(todrop, sparebits) # how many times we must apply `redrawscale` to get within reach of `mineps`
remainingbits = leftover_bottom + usebits # bits remaining beyond the lowest usable scale
f = ntuple(Returns(zero(T)), N) # the values we plan to return
good = ntuple(Returns(false), N) # track which lanes are valid
for _ in 1:numscales
scale *= redrawscale # prepare to evaluate this scale
u = ntuple(_->rand(rng, U), N) # TODO: need a SIMD version of this
willbegood = u .>= redrawbelow # values that will be accepted at this scale
fnew = trunctofloat.(u) .* scale # compute the values
f = ifelse.(good, f, fnew) # replace bad values from previous iterations
good = good .| willbegood # mark values to keep
all(good) && return f # all values are valid
end
# if we made `scale` smaller it would underflow to 0
u = ntuple(_->rand(rng, U), N) # TODO: need a SIMD version of this
u = u .>>> (8sizeof(U) - remainingbits) # we can't use all the bits at this level
fnew = trunctofloat.(u) .* mineps # compute the values
f = ifelse.(good, f, fnew) # replace bad values from previous iterations
return f
end
```
I verified this on `Float16` (that it produces all values and with appropriate frequencies), although `Float16` should maybe be done slightly differently due to the small amount of entropy and high redraw probability of 2^{-5} (which is exacerbated by blockwise computations). Hopefully I have everything set to generalize to longer floats correctly, although they require too many samples to verify emprically.
That said, this algorithm is somewhat poor in its worst-case use of entropy. It discards values with probabilities 2^{-8} for `Float32`, but when it does it utilizes only 8 bits of entropy from that draw (2^{-11} and 11 bits for `Float64`). This is probably a tradeoff, in that it streamlines the logic at the final level (which any level is likely to be) while having a larger-than-necessary worst-case number of iterations. I tried a “sig-exp”-style algorithm but it was slower and the frequencies of the very least-likely values were incorrect.
There’s likely still room for appreciable optimization, but here are some very crude benchmarks. Note that this is still missing SIMD `rand(UIntXX)` calls so I also compare it against our current scalar `rand(FloatXX)`. But how this performs with actual SIMD RNG is the actual question.
```julia-repl
julia> using BenchmarkTools
julia> @btime [block_randf(Float32,Val(16)) for _ in 1:2^14];
440.900 μs (2 allocations: 1.00 MiB)
julia> @btime [rand(Float32) for _ in 1:16*2^14];
294.200 μs (2 allocations: 1.00 MiB)
julia> @btime rand(Float32, 16*2^14);
61.900 μs (2 allocations: 1.00 MiB)
```
---
### Author:  [@dlakelan](https://discourse.julialang.org/u/dlakelan)
#### Post date: [November 1, 2023, 9:24pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/348 "2023-11-01T21:24:33Z")
Well, just for posterity and because it took a few minutes to make the graph, and because the graph is a good laugh…
Using:
```julia
function float32conv()
r = 2.0f0
exp = 30
while r < 0.0f0 || r >= 1.0
i = rand(UInt64)
a = (i >>> 32) % UInt32
b = reinterpret(Int32,i % UInt32)
r = (ldexp(Float32(a),-32) + (ldexp(Float32(b),-31-exp) - ldexp(1.0f0,-exp)))/(1.0f0 - ldexp(1.0f0,-exp+1))
end
r
end
```
And comparing it to ldexp(Float32(rand(UInt32)),-32) as well as rand(Float32)… in the far left tail below 1e-7
SO I’m looking forward to the perfect N=32 or better SIMD versions… keep up the great work @mikmoore
---
### Author:  [@mbauman](https://discourse.julialang.org/u/mbauman)
#### Post date: [November 1, 2023, 10:44pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/349 "2023-11-01T22:44:50Z")
OK, here’s another flavor of an equal-bit-width rounding down algorithm that seems really promising performance-wise. And it’s super simple. You do a bit of bit-twiddling so only the least-significant set bit remains (an exact power of two!), and then convert _that_ to a Float. So you can easily grab the resulting exponent! Make that negative at the right offset and we’re done.
```julia-auto
function f32bits(r::UInt32)
last_bit = r & -r
exponent = UInt32(253)<<23 - reinterpret(UInt32, Float32(last_bit))
exponent *= !iszero(r)
fraction = (r ⊻ last_bit)>>9 # unset the trailing one and shift
return reinterpret(Float32, exponent | fraction)
end
```
The only expensive part here is a single Int-\>Float conversion — which we already do! I see about the same speed as @mikmoore’s algorithms on my M1… but the simplicity here made me curious about Intel. With a cloud Intel AVX512 system, I see **no performance hit** … actually with _either_ Mikmoore’s algorithm or mine. And I don’t even need to do the LLVM code thing there (although it is helpful on my Mac for some unknown reason; and with it it’s slightly faster than the LLVM-ized version of Mikmoore’s alg).
I’ve not yet fiddled with iteratively adding more bits at the low-end, but I don’t think that should be too tough to add directly into `XoshiroSimd.xoshiro_bulk_simd`.
> **REPL sessions**
>
> #### Intel with AVX
>
> ```julia-repl
> julia> using Random, BenchmarkTools
>
> julia> v = Vector{Float32}(undef, 2^24);
>
> julia> @benchmark rand!($v)
> BenchmarkTools.Trial: 904 samples with 1 evaluation.
> Range (min … max): 5.120 ms … 8.815 ms ┊ GC (min … max): 0.00% … 0.00%
> Time (median): 5.480 ms ┊ GC (median): 0.00%
> Time (mean ± σ): 5.518 ms ± 257.506 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
>
> ▁▅▅▆▃▅▆█▇▆▃▅▇▄▇▁▃▄▆▂ ▂ ▃
> ▃▂▁▃▄▄▆▇▇████████████████████▆█▇▇██▅▇▇▇▄▆▅▃▄▃▃▃▃▃▃▂▂▃▂▁▁▁▃▂ ▅
> 5.12 ms Histogram: frequency by time 6.12 ms <
>
> Memory estimate: 0 bytes, allocs estimate: 0.
>
> julia> @eval Random.XoshiroSimd begin
> @inline _bits2float(x::UInt64, ::Type{Float32}) = _f32bits(x)
>
> @inline function _f32bits(r::UInt32)
> last_bit = r & -r
> exponent = UInt32(253)<<23 - reinterpret(UInt32, Float32(last_bit))
> exponent *= !iszero(r)
> fraction = (r ⊻ last_bit)>>9 # unset the trailing one and shift
> return exponent | fraction
> end
> @inline function _f32bits(x::UInt64)
> ui = (x>>>32) % UInt32
> li = x % UInt32
> return (UInt64(_f32bits(ui)) << 32) | UInt64(_f32bits(li))
> end
> @inline _f32bits(x::VecElement{UInt64}) = _f32bits(x.value)
> end
> _f32bits (generic function with 3 methods)
>
> julia> for N in [4,8,16]
> VT = :(NTuple{$N, VecElement{UInt64}})
> @eval Random.XoshiroSimd @inline _bits2float(x::$VT, ::Type{Float32}) = map(_f32bits, x)
> end
>
> julia> @benchmark rand!($v)
> BenchmarkTools.Trial: 946 samples with 1 evaluation.
> Range (min … max): 5.048 ms … 13.644 ms ┊ GC (min … max): 0.00% … 0.00%
> Time (median): 5.169 ms ┊ GC (median): 0.00%
> Time (mean ± σ): 5.275 ms ± 499.710 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
>
> ▆███▇▆▅▄▃▃▂▂▁▁
> ██████████████▇█▇▆▆▅▅▆▆▁▁▆▆▆▄▄▆▆▆▁▁▄▅▆▁▁▅▁▁▅▅▄▁▁▄▁▁▁▁▁▁▁▁▁▄ █
> 5.05 ms Histogram: log(frequency) by time 6.89 ms <
>
> Memory estimate: 0 bytes, allocs estimate: 0.
>
> ```
>
> #### Mac M1
>
> ```julia-repl
> julia> using Random, BenchmarkTools
>
> julia> v = Vector{Float32}(undef, 2^24);
>
> julia> @benchmark rand!($v)
> BenchmarkTools.Trial: 942 samples with 1 evaluation.
> Range (min … max): 5.225 ms … 6.515 ms ┊ GC (min … max): 0.00% … 0.00%
> Time (median): 5.275 ms ┊ GC (median): 0.00%
> Time (mean ± σ): 5.311 ms ± 125.292 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
>
> ▂█▇▆▆▅▅▅▄▂▁▁▁ ▁
> ██████████████▇▅█▇▆▇██▅▅▆▅▅▆▆▅▆▅▁▁▆▄▄▁▄▁▄▁▁▄▁▁▁▅▄▁▄▁▁▁▁▁▁▁▄ █
> 5.23 ms Histogram: log(frequency) by time 5.93 ms <
>
> Memory estimate: 0 bytes, allocs estimate: 0.
>
> julia> @eval Random.XoshiroSimd begin
> @inline _bits2float(x::UInt64, ::Type{Float32}) = _f32bits(x)
>
> @inline function _f32bits(r::UInt32)
> last_bit = r & -r
> exponent = UInt32(253)<<23 - reinterpret(UInt32, Float32(last_bit))
> exponent *= !iszero(r)
> fraction = (r ⊻ last_bit)>>9 # unset the trailing one and shift
> return exponent | fraction
> end
> @inline function _f32bits(x::UInt64)
> ui = (x>>>32) % UInt32
> li = x % UInt32
> return (UInt64(_f32bits(ui)) << 32) | UInt64(_f32bits(li))
> end
> @inline _f32bits(x::VecElement{UInt64}) = _f32bits(x.value)
> end
> _f32bits (generic function with 3 methods)
>
> julia> for N in [4,8,16]
> VT = :(NTuple{$N, VecElement{UInt64}})
> @eval Random.XoshiroSimd @inline _bits2float(x::$VT, ::Type{Float32}) = map(_f32bits, x)
> end
>
> julia> @benchmark rand!($v)
> BenchmarkTools.Trial: 569 samples with 1 evaluation.
> Range (min … max): 8.646 ms … 13.928 ms ┊ GC (min … max): 0.00% … 0.00%
> Time (median): 8.689 ms ┊ GC (median): 0.00%
> Time (mean ± σ): 8.792 ms ± 399.898 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
>
> █▇▅▅▃▁▁
> ████████▇▇▇█▆▆▆▇▅▄▆▄▅▅▁▆▅▅▄▁▁▄▄▁▁▄▄▁▁▁▁▁▄▁▁▁▄▁▁▁▁▁▄▅▅▁▁▄▁▁▄ ▇
> 8.65 ms Histogram: log(frequency) by time 10.2 ms <
>
> Memory estimate: 0 bytes, allocs estimate: 0.
>
> julia> for N in [4,8,16]
> VT = :(NTuple{$N, VecElement{UInt64}})
> code = """
> %i1 = bitcast <$N x i64> %0 to <$(2N) x i32>
> %i2 = sub <$(2N) x i32> zeroinitializer, %i1
> %i3 = and <$(2N) x i32> %i1, %i2
> %i4 = uitofp <$(2N) x i32> %i3 to <$(2N) x float>
> %i5 = bitcast <$(2N) x float> %i4 to <$(2N) x i32>
> %i7 = sub <$(2N) x i32> <$(join(fill("i32 2122317824", 2N), ", "))>, %i5
> %i8 = icmp eq <$(2N) x i32> %i1, zeroinitializer
> %i9 = select <$(2N) x i1> %i8, <$(2N) x i32> zeroinitializer, <$(2N) x i32> %i7
> %i10 = xor <$(2N) x i32> %i1, %i3
> %i11 = lshr <$(2N) x i32> %i10, <$(join(fill("i32 9", 2N), ", "))>
> %i12 = or <$(2N) x i32> %i9, %i11
> %i13 = bitcast <$(2N) x i32> %i12 to <$N x i64>
> ret <$N x i64> %i13
> """
> @eval Random.XoshiroSimd @inline _bits2float(x::$VT, ::Type{Float32}) = llvmcall($code, $VT, Tuple{$VT}, x)
> end
>
> julia> @benchmark rand!($v)
> BenchmarkTools.Trial: 689 samples with 1 evaluation.
> Range (min … max): 7.152 ms … 10.568 ms ┊ GC (min … max): 0.00% … 0.00%
> Time (median): 7.201 ms ┊ GC (median): 0.00%
> Time (mean ± σ): 7.260 ms ± 205.569 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
>
> ▁██▆▄▁
> ███████▇▇▇▄▆▄▆▆▄▆▄▆▄▄▆▆▅▅▅▁▄▄▄▁▄▄▁▅▁▅▅▆▄▄▅▄▄▆▆▆▆▁▅▁▁▄▁▄▁▅▁▄ ▇
> 7.15 ms Histogram: log(frequency) by time 8.01 ms <
>
> Memory estimate: 0 bytes, allocs estimate: 0.
>
> ```
---
### Author:  [@Dan](https://discourse.julialang.org/u/Dan)
#### Post date: [November 1, 2023, 11:40pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/350 "2023-11-01T23:40:03Z")
> [@mbauman](#):
>
> ```julia
> function f32bits(r::UInt32)
> last_bit = r & -r
> exponent = UInt32(253)<<23 - reinterpret(UInt32, Float32(last_bit))
> exponent *= !iszero(r)
> fraction = (r ⊻ last_bit)>>9 # unset the trailing one and shift
> return reinterpret(Float32, exponent | fraction)
> end
>
> ```
Is `f32bits` much different from:
```julia
function f32bitsB(r::UInt32)
exponent = UInt32((UInt32(126)-trailing_zeros(r)) << 23)
return reinterpret(Float32, (r >> 9)|exponent)
end
```
Perhaps `trailing_zeros` which is equivalent to conversion of a single bit into Float32 (but done in hardware pre-AVX512) is indeed better with the conversion.
But a differnet issue may be the reusing of random bits in the case `r` has more than 9 trailing bits, in which case it is both in `fraction` (regardless of xoring) and influences the exponent?
---
### Author:  [@dlakelan](https://discourse.julialang.org/u/dlakelan)
#### Post date: [November 1, 2023, 11:42pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/351 "2023-11-01T23:42:01Z")
> [@Dan](#):
>
> But a differnet issue may be the reusing of random bits in the case `r` has more than 9 trailing bits, in which case it is both in `fraction` (regardless of xoring) and influences the exponent?
Yes this was what I was trying to figure out whether there was something going on like this.
In the general case if you want to count zeros to get an exponential random variable you’ll need up to around 32 of them, and THEN the 24 bit mantissa, so I think you can do the same thing with a UInt64 broken into 2 pieces but I don’t think you can do otherwise.
---
### Author:  [@mbauman](https://discourse.julialang.org/u/mbauman)
#### Post date: [November 1, 2023, 11:53pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/352 "2023-11-01T23:53:28Z")
Those extra trailing zeros just end up as trailing zeros in the mantissa — for very small exponents. In other words, this represents the region where not all floating point values are represented. The math works out such that the steps between values in this region are constant. In the limit, there’s only an all-zero mantissa with an exponent of -32, one (most significant) bit with -31, two bits with -30, etc. You can see that the steps between all those values is 2^-32.
I had fretted about this upthread myself, but I even collected all values of a trailing zeros version to prove that it’s equivalent to the exact round-down division. I didn’t do the same test here yet.
---
### Author:  [@dlakelan](https://discourse.julialang.org/u/dlakelan)
#### Post date: [November 2, 2023, 12:11am UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/353 "2023-11-02T00:11:20Z")
> [@mbauman](#):
>
> for very small exponents.
Isn’t it all exponents less than -9 ? (I guess your point is -9-24 = -33 ?)
---
### Author:  [@danielwe](https://discourse.julialang.org/u/danielwe)
#### Post date: [November 2, 2023, 1:06am UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/354 "2023-11-02T01:06:10Z")
To make sure I understand: by using the last set bit rather than the first set bit to determine the exponent, the number you’re effectively dividing by 2^{32} and rounding here is not `r`, but rather the integer you would obtain by chopping off the last set bit and trailing zeros, and prepending one set bit and an equal number of leading zeros? Or, in case bitshifts are easier to read than words,
```julia
u = ((UInt32(1) << 31) | (r >> 1)) >> trailing_zeros(r)
```
---
### Author:  [@danielwe](https://discourse.julialang.org/u/danielwe)
#### Post date: [November 2, 2023, 1:27am UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/355 "2023-11-02T01:27:08Z")
You can use this to compare against the simple and slow conversion and confirm that your algorithm does what it intends to do:
```julia-repl
julia> slowf32bits(u::UInt32) = Float32(u, RoundDown) * 2f0^(-32);
julia> function samef32bits(r::UInt32)
u = ((UInt32(1) << 31) | (r >> 1)) >> trailing_zeros(r)
return f32bits(r) === slowf32bits(u)
end;
julia> all(samef32bits, typemin(UInt32):typemax(UInt32))
true
```
---
### Author:  [@danielwe](https://discourse.julialang.org/u/danielwe)
#### Post date: [November 2, 2023, 2:18am UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/356 "2023-11-02T02:18:04Z")
Is it possible, via some llvm magic or otherwise, to directly set the rounding mode for a single instruction? As far as I can tell, what we’re trying to do here actually has hardware support, it’s just that the global rounding mode set by julia isn’t what we need. Specifically, I can be a bit hacky and do the following to beat every other suggestion for single UInt32 → Float32:
```julia-repl
julia> Base.Rounding.setrounding_raw(Float32, Base.Rounding.to_fenv(RoundDown))
0
julia> naivef32bits(u::UInt32) = Float32(u) * 2f0^(-32)
naivef32bits (generic function with 1 method)
julia> @btime naivef32bits(rand(UInt32));
3.182 ns (0 allocations: 0 bytes)
julia> @btime f32bits(rand(UInt32)); # @mbauman's solution
4.037 ns (0 allocations: 0 bytes)
julia> @btime rand(Float32); # Status quo
2.942 ns (0 allocations: 0 bytes)
```
Changing the global rounding mode inside the conversion function is no solution, it would be both slow and not threadsafe, but if you could somehow make just this function invoke the round-down version of the hardware instruction (seems to be `vcvtsi2ss` for x86)? I’m way out of my depth here, but would love to learn.
(Didn’t include the code here to reduce clutter, but I did confirm that `naivef32bits` rounds as desired for every UInt32 after the `setrounding_raw` hack.)
---
### Author:  [@dlakelan](https://discourse.julialang.org/u/dlakelan)
#### Post date: [November 2, 2023, 4:36am UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/357 "2023-11-02T04:36:37Z")
So, just for giggles I wrote a simple function to generate 8 floats at a time, in hopes that Julia would SIMD the heck out of it automatically. It generates 4 UInt64 and chunks out 32 bits for each float, then ldexp.(mants,-32)
```julia
function eightfloat32s()
mask32 = 0xffffffff
exponents = -32
mants = (rand(UInt64),rand(UInt64),rand(UInt64),rand(UInt64))
mantissas = (mants[1] & mask32, (mants[1] >> 32) & mask32,
mants[2] & mask32, (mants[2] >> 32) & mask32,
mants[3] & mask32, (mants[3] >> 32) & mask32,
mants[4] & mask32, (mants[4] >> 32) & mask32)
ldexp.(Float32.(mantissas),exponents)
end
julia> @btime(eightfloat32s())
13.673 ns (0 allocations: 0 bytes)
(0.34529638f0, 0.38367286f0, 0.97725815f0, 0.40670988f0, 0.91882336f0, 0.12232954f0, 0.6455661f0, 0.87967235f0)
13.673/8
1.709125
```
So that’s 1.7 nanoseconds per Float32 with relative simple code.
```julia
(let a = zeros(Float32,8)
@btime (rand!($a,Float32))
end)
12.295 ns (0 allocations: 0 bytes)
julia> 13.673/12.295
1.1120780805205368
```
So about an 11 percent hit compared to rand! for 8 at a time…
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [November 2, 2023, 8:28am UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/358 "2023-11-02T08:28:44Z")
> [@dlakelan](#):
>
> It comes down to a question of frequency in real world code.
I disagree with this. If something is a bug, it should be fixed, regardless of the frequency. If it isn’t, we just document the API to inform users.
Practically, I think that “may output 0” is easier to live with than “cannot output 0”, as it provides more flexibility for implementations, which may or may not use it.
But I think that 0 is a side issue in this topic. And I am no longer sure what the real issue is: for `rand(Float64)`, basically anything goes, including inverse transform sampling, because it is dense enough for practical purposes as it is, and IMO the current implementation is fine.
For `rand(Float32)` this is no longer true, but that’s mainly because of `Float32`, not `rand`. Fixing problems that arise from nonlinear mappings by changing `rand` seems very roundabout to me. I think we should just have the current expectations of random number generators documented.
---
### Author:  [@Dan](https://discourse.julialang.org/u/Dan)
#### Post date: [November 2, 2023, 9:07am UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/359 "2023-11-02T09:07:01Z")
> [@Tamas\_Papp](#):
>
> If something is a bug, it should be fixed
This is agreed.
I could post a poll: _Should `rand` return every Float32 in the range?_ but there is no need. The result is quite predictable: 80% or more will say no. This is driven a lot by _status quo bias_ which is strong in the community (and a topic for a different thread). Given this poll conclusion, and that currently there are numbers in the range and not outputted by `rand`, any such number has _at least_ the same claim to appear as `0`.
If we have the equivalent of `ifelse(res == 0, nextfloat(res), res)` it would cost little in performance, and would introduce a subnormal which all the non-full-range-output status quo algos would not produce anyway.
As for:
> [@dlakelan](#):
>
> It comes down to a question of frequency in real world code
Real world code _is_ relevant as we don’t have a specification for `rand` and thus need informative signals wherever we find them.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [November 2, 2023, 9:25am UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/360 "2023-11-02T09:25:31Z")
> [@Dan](#):
>
> it would cost little in performance, and would introduce a subnormal
I still do not see the motivation for this though. Please understand that I can conceive multiple scenarios where I would be bitten by non-linear transformations of a `rand(Float32)`, but in all of them I would just do something else in the first place, eg use `rand(Float64)` in general, but if I really have to use `rand(Float32)` and I am interested in the tail (eg extreme values), just scale down draws, eg by 2^8 or 2^{16}, depending on what I am interested in, and transform from those.
> [@Dan](#):
>
> Real world code _is_ relevant
Real world _problems_ would be more interesting for me. Specifically, a description of _what_ the user is trying to achieve, not _how_.
I think changes to the language and the standard libraries should ideally be motivated by these. Especially because if we make something a requirement of `rand` _in general_, all implementations should support it. This is not a trivial change, changing the code in `Random` is only the start of it.
---
### Author:  [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)
#### Post date: [November 2, 2023, 11:00am UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/361 "2023-11-02T11:00:21Z")
> [@danielwe](#):
>
> Is it possible, via some llvm magic or otherwise, to directly set the rounding mode for a single instruction?
Not on x86.
The rounding mode is controlled by some control register. We can change the rounding mode, either via machine-dependent assembly (bad idea) or via some `ccall`, e.g. to libc. Changing that is unlike other register changes (e.g. `movq`) that are somewhere between cheap and free – it’s a rare probably micro-coded event for the CPU.
We can do that for the SIMD code, since this amortizes over a large array. We cannot do that for a single `rand()` call.
> [@Tamas\_Papp](#):
>
> I still do not see the motivation for this though. Please understand that I can conceive multiple scenarios where I would be bitten by non-linear transformations of a `rand(Float32)`, but in all of them I would just do something else in the first place, eg use `rand(Float64)` in general, but if I really have to use `rand(Float32)` and I am interested in the tail (eg extreme values), just scale down draws, eg by 2^8282^8 or 2^{16}2162^{16}, depending on what I am interested in, and transform from those.
It hurts me to re-open that thing, but I still think about the `rand()^-2` code for sampling from pareto.
That’s not what Distributions.jl uses, but after reviewing your recommendation upthread (the James Gentle book), I feel pretty justified in my claim that this mistake does not prove that the user is hopelessly numerically illiterate and that we can throw them under the bus (which was your position, I think?).
I mean, the James Gentle book
1. recommends to use that algorithm (i.e. inverse cfd on `U(0,1)`) in the section about Pareto. I believe that James Gentle is very aware of what non-linear transforms with singularities can be done on floating point numbers without blowing up.
2. acknowledges in the introduction that sampling from `U(0,1)` cannot be adequately done the way we do it (e.g. by sampling uniform from (1,2) and subtracting 1) due to all the issues with loss of precision close to zero we talked about, and
3. recommends in the exercises and examples to use standard library code for `U(0,1)` (apparently either without having read said standard library code, or without noticing the contradiction to point 2)
Can we really throw people under the bus who follow that book? Or is the book just non-representative of the scientific literature, and we 100% should tell people to read something better? Do you want to retract your recommendation and recommend another book? (preferably textbook-level, not monograph)
I’m not asking you for the recommendation because I’m too lazy to google something, or to score a gotcha point.
It’s because I lack the necessary subtext and community and implied qualificactions (it’s not my field! And I’m out of the academy anyways) to adequately grab something. A la “hey, you must read preprint xyz! It is wrong and will be retracted soon, but the idea so so cool and it’s the most credible attack on conjecture xyz in years”, or “take a look at book xyz; especially the chapter on abc is good. The parts about xyz are really bad though, go somewhere else for that”.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [November 2, 2023, 11:17am UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/362 "2023-11-02T11:17:30Z")
> [@foobar\_lv2](#):
>
> I mean, the James Gentle book
>
> 1. recommends to use that algorithm (i.e. inverse cfd on `U(0,1)`) in the section about Pareto. I believe that James Gentle is very aware of what non-linear transforms with singularities can be done on floating point numbers without blowing up
The implicit context for this is using `Float64`. If you do that, the inverse cdf transform works just fine for practical purposes.
> [@foobar\_lv2](#):
>
> Can we really throw people under the bus who follow that book?
The book was written in 2003 (if I remember correctly). `Float64` had been well-supported at that point on the CPU level for more than a decade, so it has become the de facto standard for numerical work — most textbooks use it implicitly, without even mentioning 32-bit floats.
Working with 32-bit floats is wrought with a lot of complications, not restricted to random numbers. If the user insists on 32-bit floats, I would say that the burden is on them to become informed about its consequences.
Specifically in Julia, you cannot run into 32-bit floats by accident, you have to request them explicitly and then make sure that you remain in that precision with carefully written generic code (eg `eps(T)` instead of `eps()`, etc).
(Can we please dispense with the “innocents suffer because `rand(Float32)` is not dense enough around zero” rhetoric? It is really not moving this discussion forward.)
> [@foobar\_lv2](#):
>
> Do you want to retract your recommendation and recommend another book?
No, it is a fine textbook. The problem is that you insist on `Float32` computations. Those were a thing in the early days of computing, you will find a ton of papers if you have a specific problem to solve (which, incidentally, I am still missing from this discussion).
[Previous page](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184.md?page=17)
[Next page](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184.md?page=19)