Correct bias in random float sampling.

Fixes JuliaLang#33222.

Also changed some hardcoded test values accordingly.

Author: tpapp

stdlib/Random/src/Random.jl

@@ -309,8 +309,8 @@ Pick a random element or array of random elements from the set of values specifi
 * an `AbstractDict` or `AbstractSet` object,
 * a string (considered as a collection of characters), or
 * a type: the set of values to pick from is then equivalent to `typemin(S):typemax(S)` for
-  integers (this is not applicable to [`BigInt`](@ref)), to ``[0, 1)`` for floating
-  point numbers and to ``[0, 1)+i[0, 1)]`` for complex floating point numbers;
+  integers (this is not applicable to [`BigInt`](@ref)), to ``(0, 1)`` for floating
+  point numbers and to ``(0, 1)+i(0, 1)]`` for complex floating point numbers;
 
 `S` defaults to [`Float64`](@ref).
 

stdlib/Random/src/generation.jl

@@ -22,17 +22,26 @@ Sampler(::Type{RNG}, ::Type{T}, n::Repetition) where {RNG<:AbstractRNG,T<:Abstra
 # generic random generation function which can be used by RNG implementors
 # it is not defined as a fallback rand method as this could create ambiguities
 
-rand(r::AbstractRNG, ::SamplerTrivial{CloseOpen01{Float16}}) =
-    Float16(reinterpret(Float32,
-                        (rand(r, UInt10(UInt32)) << 13)  | 0x3f800000) - 1)
+function rand(r::AbstractRNG, ::SamplerTrivial{CloseOpen01{Float16}})
+    # bias correction below is equivalent to adding eps(Float16)/2, see #33222
+    z = reinterpret(Float32, (rand(r, UInt10(UInt32)) << 13)  | 0x3f800000)
+    Float16(z - prevfloat(Float16(1)))
+end
 
-rand(r::AbstractRNG, ::SamplerTrivial{CloseOpen01{Float32}}) =
-    reinterpret(Float32, rand(r, UInt23()) | 0x3f800000) - 1
+function rand(r::AbstractRNG, ::SamplerTrivial{CloseOpen01{Float32}})
+    # bias correction below is equivalent to adding eps(Float32)/2, see #33222
+    reinterpret(Float32, rand(r, UInt23()) | 0x3f800000) - prevfloat(1f0)
+end
 
-rand(r::AbstractRNG, ::SamplerTrivial{CloseOpen12_64}) =
+function rand(r::AbstractRNG, ::SamplerTrivial{CloseOpen12_64})
+    # NOTE this has a small bias of eps()/2, see discussion in #33222
     reinterpret(Float64, 0x3ff0000000000000 | rand(r, UInt52()))
+end
 
-rand(r::AbstractRNG, ::SamplerTrivial{CloseOpen01_64}) = rand(r, CloseOpen12()) - 1.0
+function rand(r::AbstractRNG, ::SamplerTrivial{CloseOpen01_64})
+    # bias correction below is equivalent to adding eps()/2, see #33222
+    rand(r, CloseOpen12()) - prevfloat(1.0)
+end
 
 #### BigFloat
 

stdlib/Random/test/runtests.jl

@@ -33,10 +33,10 @@ end
 @test length(randn(ComplexF64, 4, 5)) == 20
 @test length(bitrand(4, 5)) == 20
 
-@test rand(MersenneTwister(0)) == 0.8236475079774124
-@test rand(MersenneTwister(42)) == 0.5331830160438613
+@test rand(MersenneTwister(0)) == 0.8236475079774125
+@test rand(MersenneTwister(42)) == 0.5331830160438614
 # Try a seed larger than 2^32
-@test rand(MersenneTwister(5294967296)) == 0.3498809918210497
+@test rand(MersenneTwister(5294967296)) == 0.3498809918210498
 
 # Test array filling, Issues #7643, #8360
 @test rand(MersenneTwister(0), 1) == [0.8236475079774124]
@@ -307,7 +307,8 @@ let mt = MersenneTwister(0)
         Random.seed!(mt, 0)
         rand(mt) # this is to fill mt.vals, cf. #9040
         rand!(mt, A) # must not segfault even if Int(pointer(A)) % 16 != 0
-        @test A[end-4:end] == [0.3371041633752143, 0.41147647589610803, 0.6063082992397912, 0.9103565379264364, 0.16456579813368521]
+        @test A[end-4:end] == [0.3371041633752143, 0.41147647589610803, 0.6063082992397912,
+                               0.9103565379264366, 0.16456579813368533]
     end
 end