Proposal: randfloat() for sampling from _all_ floats in [0,1)?

[milankl]

As it's written in the documentation, the standard approach of for rand(Float32)/rand(Float64)

An UInt64 will firstly be converted to a Float64 that is perfectly uniformly distributed in [1.0, 2.0), and then minus one. This is a very fast approach, but not completely ideal, since the one bit is wasted. The current default RNG in Base.Random library does the same thing, so it also causes some tricky problems

is probably good enough in most cases, but not perfect. For Float32 it samples only from 2^23 = 8,388,608 floats in [+0,1) (as that’s the number of floats in [1,2)) whereas there are 1,065,353,216 (=2^30-2^23) in that range in total. So rand(Float32) samples from about 1% possible floats, and rand(Float64) from ~0.1% (2^52/(2^62-2^52)). I propose a function that solves that problem by drawing the exponent bits at correct chances first and then combine them with significant bits, similar to the following

using RandomNumbers.Xorshifts
const Xoroshiro128P = Xoroshiro128Plus()

"""Random number generator for floats in [0,1) that samples from *all* floats.""" 
function randfloat(::Type{Float32})
    # create exponent bits in 0000_0000 to 0111_1111
    # at following chances
    # e=01111110 at 50.0% for [0.5,1.0)
    # e=01111101 at 25.0% for [0.25,0.5)
    # e=01111100 at 12.5% for [0.125,0.25)
    # ... etc
    ui = rand(Xoroshiro128P,UInt128) >> 1 | 0x1 # first bit always 0, last always 1
    lz = leading_zeros(ui) % UInt8              # 2 leading zeros is 50%, 3 is 25%, etc.  
    e = ((0x7f - lz) % UInt32) << 23            # convert lz to exponent bits of float32

    # significant bits
    f = rand(Xoroshiro128P,UInt32) >> 9 
    
    # combine exponent and significand (sign always 0)
    return reinterpret(Float32,e | f)
end

I don't know whether something like this exists already in Julia somewhere, I don't even know whether this is an established algorithm. I thought it would be nice to have this function and different to implemented here, it could take the actual RNG as an argument. I've just used Xoroshiro128+ as it's fast. I'm happy to implement a Float64 version too, this requires technically a 1024bit UInt, but in practice one could first draw a UInt64 and only if 0x0 is drawn (does that actually happen?) one draws another UInt64 and count the leading zeros again. The above implementation seems to be reasonably fast, about 3x slower than Xoroshiro128+ with the "[1,2)-1"-technique but still faster than Julia's Base rand(Float32):

julia> @btime randfloat($Float32)   # this approach
  4.603 ns (0 allocations: 0 bytes)
julia> @btime rand($Float32)          # Julia's Base MT19337
  6.675 ns (0 allocations: 0 bytes)
julia> @btime rand($Xoroshiro128P,$Float32)    # Xoroshiro128+ with [1,2)-1 technique
  1.544 ns (0 allocations: 0 bytes)

Again, there might be ways to speed this up further by first drawing a UInt64 and only if that's 0x0 draw another one.

[milankl]

@eapax, I finally managed to write an algorithm for this problem in essentially 3 lines :)

[milankl]

Okay, rand(UInt128) introduces a bit of overhead, doing rand(UInt64) instead shaves off another 1.5ns

ui = rand(Xoroshiro128P,UInt64) >> 1    # first bit always 0
# count leading zeros of UInt128 in two steps: 2. only if 
lz = ui == 0 ?  # 2 leading zeros is at 50%, 3 at 25%, etc.
    64+leading_zeros(rand(Xoroshiro128P,UInt64) | 0x1) : leading_zeros(ui) 
e = ((127 - lz) % UInt32) << 23     # convert lz to exponent bits of float32

then yields

julia> @btime randfloat($Float32)
  3.060 ns (0 allocations: 0 bytes)

So 2x faster than MersenneTwister and only 2x slower than rand(Xoroshiro128Plus,Float32). The "I don't care about number smaller than 5e-20"-version which avoids the ?-branching is also with 2.75ns probably only irrelevantly faster.

[milankl]

Turned this into a PR https://github.com/sunoru/RandomNumbers.jl/pull/73

[sunoru]

Looks great! Thank you very much for suggesting this method to sample from all floats in [0,1).

I made some comments in your PR, and can you also try to make the CI tests pass?

[milankl]

Can you comment on whether it's possible with, say, Xoroshiro128Plus (but also other RNGs) to hit UInt64(0)? I'm not sure whether that's the reason why it's not allowed as seed.

[milankl]

After ui = rand(RNG,UInt64) >> 1 the ui==0 case can also be hit with ui=0x1, so for Float32 that should be fine. For Float64, however, the idea is to sample several ui = rand(RNG,UInt64) and only sample the next when the previous one is ui==0 but would that ever happen?

[milankl]

That's also a nice quality test

julia> using RandomNumbers
julia> x = rand(Float32,100_000_000);
julia> y = randfloat(Float32,100_000_000);
julia> length(Set(x))
8388552

julia> length(Set(y))
32781328

2^23=8388608 is the size of the set that rand(Float32) draws from, so on when you draw 1e8 float32s then you hit almost all of them and on average each of them 10 times. This problem is clearly mitigated with randfloat and probably as good as possible

julia> z = Float32.(rand(100_000_000));

julia> length(Set(z))
32784392

[sunoru]

Can you comment on whether it's possible with, say, Xoroshiro128Plus (but also other RNGs) to hit UInt64(0)? I'm not sure whether that's the reason why it's not allowed as seed.

All RNGs in xorshift family cannot hit UInt64(0), since the rotation and xor operations on zero values do not have any effects. This is also the reason why it's not allowed as seed, since there is not much difference between a seed and a state in xorshift.

After ui = rand(RNG,UInt64) >> 1 the ui==0 case can also be hit with ui=0x1, so for Float32 that should be fine. For Float64, however, the idea is to sample several ui = rand(RNG,UInt64) and only sample the next when the previous one is ui==0 but would that ever happen?

I am not sure if we really need to resample a random number if we hit ui==0. The probability is just too small. As a result, a large number of numbers between 0-2^-64, around 2^23 * 63 = 528482304 states (for Float32) have to be wasted but I think it is OK?

[sunoru]

And the quality test looks good!

[milankl]

I am not sure if we really need to resample a random number if we hit ui==0. The probability is just too small. As a result, a large number of numbers between 0-2^-64, around 2^23 * 63 = 528482304 states (for Float32) have to be wasted but I think it is OK?

Yes, but as I said earlier

The "I don't care about number smaller than 5e-20"-version which avoids the ?-branching is also with 2.75ns probably only irrelevantly faster.

However, we could also use the last 23 bits of ui for the significand, if we already say that's too unlikely to be important, which might speed up things further?

[milankl]

function randfloatx(rng::AbstractRNG,::Type{Float32})
    ui = rand(rng,UInt64) >> 1    # first bit always 0
    e = ((127 - leading_zeros(ui)) % UInt32) << 23 
    return reinterpret(Float32,e | ((ui % UInt32) & 0x007f_ffff))
end

times as

julia> @btime randfloatx($Float32)
  1.904 ns (0 allocations: 0 bytes)

That's with Xoroshiro128Plus and reusing the last 23bits for the significand. It still samples from 5e-20 to prevfloat(1f0) and only 5e-20 to 5e-13 will be affected with the reusing of the last bits for the significand.

[sunoru]

It looks correct. I am not sure if this influence could be important in some cases, but I guess it's not. And the speed is nice, maybe we can rerun the TestU01 big crush tests on it to see the statistical performance.

[milankl]

I guess that if rand(Float32)/rand(Float64) already passes TestU01 then randfloat should so too. Unless, currently there's 0 chance that randfloat produces 0f0/0. Is that important?

[milankl]

I think a better way to do things is

function randfloat(rng::AbstractRNG,::Type{Float32})
    ui = rand(rng,UInt64)
    lz = leading_zeros(ui)
    e = ((126 - lz) % UInt32) << 23
    ui = lz > 40 ? rand(rng,UInt64) : ui
    return reinterpret(Float32,e | ((ui % UInt32) & 0x007f_ffff))
end

The ?-branch doesn't introduce any overhead for me and avoids the bits to be double used. Then randfloat really samples all 64*2^23=536870912 floats between 2.7105054f-20 and prevfloat(1f0). That's 50% of all float32s in [0,1). rand(Float32), instead, only samples from 2^32=8388608 (<1%).

Similar for Float64: randfloat samples from 64x as many float64s in [0,1) than rand. That's an improvement from 0.1% to 6% of all floa64s in [0,1). Timings are now

julia> @btime rand($Xoroshiro128P,$Float64)
  1.849 ns (0 allocations: 0 bytes)

julia> @btime rand($Xoroshiro128P,$Float32)
  1.544 ns (0 allocations: 0 bytes)

with rand, vs

julia> @btime randfloat($Xoroshiro128P,$Float64)
  2.149 ns (0 allocations: 0 bytes)

julia> @btime randfloat($Xoroshiro128P,$Float32)
  1.763 ns (0 allocations: 0 bytes)

with randfloat.