Conversion to Float

[simonbyrne]

The "obvious" approach of converting to a float, then multiplying by a scale factor is slow, and has the potential to give an answer of 1.0.

The easiest option is:

import Base: significand_mask, exponent_one
f1(u::UInt64) = reinterpret(Float64, exponent_one(Float64) | significand_mask(Float64) & u) - 1.0

unfortunately this has the downside that the last bit will always be zero, so you only get 52 bits of randomness per float: see JuliaLang/julia#16344.

A slightly more advanced option (based on this proposal) is:

import Base: significand_mask, significand_bits, exponent_half
function f2(u::UInt64)
    f = reinterpret(Float64, exponent_half(Float64) | 0x001f_ffff_ffff_ffff & u)
    if (u >> significand_bits(Float64)) &1 == 1
        f-= 0.5
    end
    f
end

This gets us 53 bits, but introduces a branch. There might be some clever stuff we can do here though.

[sunoru]

I had believed what I learned in art4711/random-double#2 could be the best approach to convert UInts to Floats, but the expense of time seemed too large. See https://github.com/sunoru/notes/blob/master/GSoC%202016/Conversion%20to%20Float.ipynb for my performance test.

So I'd like to use the first (simplest and very fast) approach currently.

[chethega]

An alternative is

julia> function f3(x::UInt64)
       tz = trailing_zeros(x)
       res = reinterpret(Float64, (x>>12) | (1022-tz)<<52)
       end

This has very small bias: there is a 1/4000 chance that the least significant bit is correlated with the exponent (similar for next bits), and there is a cutoff that lumps all numbers smaller than f3(UInt(0)) = 2.710505431213761e-20 into the same bin. Otherwise, the significant bits are uniform random and the exponent is exponentially distributed, which is exactly what we want: Sample a uniform random number of infinite precision and round it to floating point. This is a very acceptable deviation and avoids all subnormals. This algorithm is faster than random generation, but does not simd, which makes it slower than f1 and f2. Nevertheless, I think this is what one should use when natively generating random integers (mersenne twister only generates 52 bit of entropy for a float, so probably should not use this).

Note that a correct PRNG with Float64 output should be able to generate 0.0 with a probability of nextfloat(0.0)=5.0e-324 which is "never ever ever", compared to the f3 truncation problem, occuring with about 1e-20 probability, i.e. "practically never".

Much more important is the conversion to Float32:

julia> function genzero()
       n=0
       while rand(Float32)!=0
       n+=1
       end
       return n
       end
julia> for i=1:30
       @show genzero()
       end
genzero() = 10667263
genzero() = 9663774
genzero() = 1169701
genzero() = 6097731
...

A correct RNG would output 0.0f with a probability of at most nextfloat(0.0f0) = 1.0f-45, which is approximately never ever, instead of every couple million draws.

[milankl]

@chethega I think your f3 and what I propose here: https://github.com/sunoru/RandomNumbers.jl/issues/72 goes in a similar direction!