Output distribution of rand(Float32) and rand(Float64), thread 2

Internals & Design
,
77

Ok, I had some time to think about the problem and came up with this, though I don’t think it’s really “better” than Float32(rand(Float64))

function rand32float()
   c = -1.0
   while c <0.0 || c >= 1.0
     i = rand(UInt64)
     a = Float32(i >> 32)/2^(32)
     b = Float32(mod(i,2^32)/2^32)
     c = (a + b/2^(12)-1/2^(12))/(1 - 1/2^(12))
    end
    return c
end

EDIT: error, above I don’t preserve Float32 type, so this fixes it I think:

function rand32float()
          c = -1.0f0
          while c <0.0f0 || c >= 1.0f0
            i = rand(UInt64)
            a = Float32(i >> 32)/2^(32)
            b = Float32(mod(i,2^32))/2^32
            c = (a + b/2^(12)-1.0f0/2^(12))/(1.0f0 - 1.0f0/2^(12))
           end
           return c
       end

Simple histograms suggest that I’m right but I didn’t write a proof.

image
image815×583 23.9 KB

On my machine:

julia> @btime rand32float()
  4.969 ns (0 allocations: 0 bytes)
0.8600827395260989

julia> @btime Float32(rand(Float64))
  2.845 ns (0 allocations: 0 bytes)
0.4434461f0

EDIT: timing of the newer version:

julia> @btime rand32float()
  4.408 ns (0 allocations: 0 bytes)
0.06446942f0

julia> @btime Float32(rand(Float64))
  2.855 ns (0 allocations: 0 bytes)
0.5801686f0
1 Like