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

<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:** [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")

</div>

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)

```

---

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