Fallback `rand(Float32/64)` discards more entropy than necessary

[danielwe]

I just noticed that the current default Xoshiro RNG generates uniform floats with more entropy than other RNGs. The reason appears to be that it has it's own 24/53-bit implementation of the conversion from random bits to uniform floats at

julia/stdlib/Random/src/Xoshiro.jl

Lines 207 to 211 in 7014681

	rand(r::Union{TaskLocalRNG, Xoshiro}, ::SamplerTrivial{CloseOpen01{Float32}}) =
	Float32(rand(r, UInt32) >>> 8) * Float32(0x1.0p-24)
	rand(r::Union{TaskLocalRNG, Xoshiro}, ::SamplerTrivial{CloseOpen01_64}) =
	Float64(rand(r, UInt64) >>> 11) * 0x1.0p-53

while the generic conversion uses a bit twiddling 23/52-bit implementation at

julia/stdlib/Random/src/generation.jl

Lines 29 to 35 in 7014681

	rand(r::AbstractRNG, ::SamplerTrivial{CloseOpen01{Float32}}) =
	reinterpret(Float32, rand(r, UInt23()) | 0x3f800000) - 1
	rand(r::AbstractRNG, ::SamplerTrivial{CloseOpen12_64}) =
	reinterpret(Float64, 0x3ff0000000000000 | rand(r, UInt52()))
	rand(r::AbstractRNG, ::SamplerTrivial{CloseOpen01_64}) = rand(r, CloseOpen12()) - 1.0

The explanation for this is here:

We still use the 52-bit version for our older RNGs so that they give the same values as before, but for the new xoshiro generator we're using shift-by-11-and-multiply currently.

Originally posted by @JeffBezanson in #33222 (comment)

However, one consequence is that the 23/52-bit conversions are also used by packages that implement their own RNG and provide only the minimal bit sampling interface, relying on generic fallbacks to produce floats. This includes, e.g., the native RNGs in CUDA.jl and StableRNGs.jl, but also any RNG library developed in the future.

Not sure if this is intended or unintended? I suppose it's desirable in the case of StableRNGs, but not so much in the case of future implementations. Wouldn't it be better to have the 24/53-bit conversion be the generic fallback, and rather specialize to the 23/52-bit conversion for legacy RNGs? StableRNGs could implement its own conversion if it wants to promise cross-version stability of float sequences.

Noticed this while investigating JuliaGPU/CUDA.jl#1464

[rfourquet]

I guess the default implementation using 52 bits comes from the fact that dSFMT's implementation of MersenneTwsister generates natively these 52 bits, but it's indeed probably a good idea to change the default, as most RNGs generate natively 64 uniform bits instead of 52. Nowadays MersenneTwister in Julia can also produce "almost natively" 64 bits (i.e. the cost per bit is not much higher for 64 than for 52), so it would not be a no-go to make the change also for this RNG if the cost is deemed justifed.