From 3a036efd773e4a6daafaa8b077980cabadedc221 Mon Sep 17 00:00:00 2001 From: Milan Date: Thu, 18 Feb 2021 21:00:12 +0000 Subject: [PATCH 1/3] randfloat(Float16) --- src/randfloat.jl | 24 +++++++++++++++++++++ test/randfloat.jl | 54 +++++++++++++++++++++++++++++++++++------------ 2 files changed, 64 insertions(+), 14 deletions(-) diff --git a/src/randfloat.jl b/src/randfloat.jl index ea97a78..7830e7e 100644 --- a/src/randfloat.jl +++ b/src/randfloat.jl @@ -59,6 +59,30 @@ function randfloat(rng::Random.AbstractRNG,::Type{Float64}) return reinterpret(Float64,e | (ui & 0x000f_ffff_ffff_ffff)) end +"""Random number generator for Float16 in [0,1) that samples from +all 15360 float16s in that range.""" +function randfloat(rng::Random.AbstractRNG,::Type{Float16}) + # create exponent bits in 00000 to 01110 + # at following chances + # e=01110 at 50.0% for [0.5,1.0) + # e=01101 at 25.0% for [0.25,0.5) + # e=01100 at 12.5% for [0.125,0.25) + # ... + ui = rand(rng,UInt32) | 0x0002_0000 + # set 15th bit to 1 to have at most 14 leading zeros. + + # count leading zeros of random UInt64 in several steps + # 0 leading zeros at 50% chance + # 1 leading zero at 25% chance + # 2 leading zeros at 12.5% chance etc. + # then convert leading zeros to exponent bits of Float16 + lz = leading_zeros(ui) + e = ((14 - lz) % UInt32) << 10 + + # combine exponent and significand (sign always 0) + return reinterpret(Float16,(e | (ui & 0x0000_03ff)) % UInt16) +end + # use stdlib default RNG as a default here too randfloat(::Type{T}=Float64) where T = randfloat(GLOBAL_RNG,T) randfloat(rng::Random.AbstractRNG) = randfloat(rng,Float64) diff --git a/test/randfloat.jl b/test/randfloat.jl index f2b42f7..1a7368d 100644 --- a/test/randfloat.jl +++ b/test/randfloat.jl @@ -5,11 +5,14 @@ using Test @test typeof(randfloat()) == Float64 @test typeof(randfloat(Float32)) == Float32 @test typeof(randfloat(Float64)) == Float64 + @test typeof(randfloat(Float16)) == Float16 # arrays + @test typeof(randfloat(Float16,2)) <: Vector{Float16} @test typeof(randfloat(Float32,2)) <: Vector{Float32} @test typeof(randfloat(Float64,2)) <: Vector{Float64} + @test typeof(randfloat(Float16,2,3)) <: Matrix{Float16} @test typeof(randfloat(Float32,2,3)) <: Matrix{Float32} @test typeof(randfloat(Float64,2,3)) <: Matrix{Float64} end @@ -18,70 +21,93 @@ end import RandomNumbers.Xorshifts.Xorshift64 rng = Xorshift64() @test typeof(randfloat(rng)) == Float64 + @test typeof(randfloat(rng,Float16)) == Float16 @test typeof(randfloat(rng,Float32)) == Float32 @test typeof(randfloat(rng,Float64)) == Float64 # arrays + @test typeof(randfloat(rng,Float16,2)) <: Vector{Float16} @test typeof(randfloat(rng,Float32,2)) <: Vector{Float32} @test typeof(randfloat(rng,Float64,2)) <: Vector{Float64} + @test typeof(randfloat(rng,Float16,2,3)) <: Matrix{Float16} @test typeof(randfloat(rng,Float32,2,3)) <: Matrix{Float32} @test typeof(randfloat(rng,Float64,2,3)) <: Matrix{Float64} import RandomNumbers.MersenneTwisters.MT19937 rng = MT19937() @test typeof(randfloat(rng)) == Float64 + @test typeof(randfloat(rng,Float16)) == Float16 @test typeof(randfloat(rng,Float32)) == Float32 @test typeof(randfloat(rng,Float64)) == Float64 # arrays + @test typeof(randfloat(rng,Float16,2)) <: Vector{Float16} @test typeof(randfloat(rng,Float32,2)) <: Vector{Float32} @test typeof(randfloat(rng,Float64,2)) <: Vector{Float64} + @test typeof(randfloat(rng,Float16,2,3)) <: Matrix{Float16} @test typeof(randfloat(rng,Float32,2,3)) <: Matrix{Float32} @test typeof(randfloat(rng,Float64,2,3)) <: Matrix{Float64} end -@testset "Distribution is uniform [0,1)" begin +@testset "Distribution is uniform [0,1) Float64" begin N = 100_000_000 - # Float64 x = randfloat(N) @test maximum(x) > 0.999999 @test minimum(x) < 1e-7 - exponent_mask = 0x7ff0_0000_0000_0000 - exponent_bias = 1023 - # histogram for exponents H = fill(0,64) for xi in x # 0.5 is in H[1], 0.25 is in H[2] etc - e = reinterpret(UInt64,xi) & exponent_mask - H[exponent_bias-Int(e >> 52)] += 1 + e = reinterpret(UInt64,xi) & Base.exponent_mask(Float64) + H[Base.exponent_bias(Float64)-Int(e >> 52)] += 1 end # test that the most frequent exponents occur at 50%, 25%, 12.5% etc. for i in 1:10 @test isapprox(H[i]/N,2.0^-i,atol=5e-4) end - - # FLOAT32 +end + +@testset "Distribution is uniform [0,1) Float32" begin + N = 100_000_000 + x = randfloat(Float32,N) @test maximum(x) > prevfloat(1f0,5) @test minimum(x) < 1e-7 - exponent_mask = 0x7f80_0000 - exponent_bias = 127 - # histogram for exponents H = fill(0,64) for xi in x # 0.5f0 is in H[1], 0.25f0 is in H[2] etc - e = reinterpret(UInt32,xi) & exponent_mask - H[exponent_bias-Int(e >> 23)] += 1 + e = reinterpret(UInt32,xi) & Base.exponent_mask(Float32) + H[Base.exponent_bias(Float32)-Int(e >> 23)] += 1 + end + + # test that the most frequent exponents occur at 50%, 25%, 12.5% etc. + for i in 1:10 + @test isapprox(H[i]/N,2.0^-i,atol=5e-4) + end +end + +@testset "Distribution is uniform [0,1) Float16" begin + N = 100_000_000 + + x = randfloat(Float16,N) + @test maximum(x) == prevfloat(one(Float16)) + @test minimum(x) == zero(Float16) # zero can be produced + + # histogram for exponents + H = fill(0,32) + for xi in x # 0.5f0 is in H[1], 0.25f0 is in H[2] etc + e = reinterpret(UInt16,xi) & Base.exponent_mask(Float16) + H[Base.exponent_bias(Float16)-Int(e >> 10)] += 1 end # test that the most frequent exponents occur at 50%, 25%, 12.5% etc. for i in 1:10 @test isapprox(H[i]/N,2.0^-i,atol=5e-4) end + println(H) end \ No newline at end of file From 52686cc8d9849c6b8ce8374fefac1cef61d9ee95 Mon Sep 17 00:00:00 2001 From: Milan Date: Thu, 18 Feb 2021 21:16:23 +0000 Subject: [PATCH 2/3] Base.exponent_mask/bias fails on 1.3 --- test/randfloat.jl | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/test/randfloat.jl b/test/randfloat.jl index 1a7368d..237c0ef 100644 --- a/test/randfloat.jl +++ b/test/randfloat.jl @@ -58,11 +58,14 @@ end @test maximum(x) > 0.999999 @test minimum(x) < 1e-7 + exponent_mask = 0x7ff0000000000000 + exponent_bias = 1023 + # histogram for exponents H = fill(0,64) for xi in x # 0.5 is in H[1], 0.25 is in H[2] etc - e = reinterpret(UInt64,xi) & Base.exponent_mask(Float64) - H[Base.exponent_bias(Float64)-Int(e >> 52)] += 1 + e = reinterpret(UInt64,xi) & exponent_mask + H[exponent_bias-Int(e >> 52)] += 1 end # test that the most frequent exponents occur at 50%, 25%, 12.5% etc. @@ -78,11 +81,14 @@ end @test maximum(x) > prevfloat(1f0,5) @test minimum(x) < 1e-7 + exponent_mask = 0x7f800000 + exponent_bias = 127 + # histogram for exponents H = fill(0,64) for xi in x # 0.5f0 is in H[1], 0.25f0 is in H[2] etc - e = reinterpret(UInt32,xi) & Base.exponent_mask(Float32) - H[Base.exponent_bias(Float32)-Int(e >> 23)] += 1 + e = reinterpret(UInt32,xi) & exponent_mask + H[exponent_bias-Int(e >> 23)] += 1 end # test that the most frequent exponents occur at 50%, 25%, 12.5% etc. @@ -98,11 +104,14 @@ end @test maximum(x) == prevfloat(one(Float16)) @test minimum(x) == zero(Float16) # zero can be produced + exponent_mask = 0x7c00 + exponent_bias = 15 + # histogram for exponents - H = fill(0,32) + H = fill(0,16) for xi in x # 0.5f0 is in H[1], 0.25f0 is in H[2] etc - e = reinterpret(UInt16,xi) & Base.exponent_mask(Float16) - H[Base.exponent_bias(Float16)-Int(e >> 10)] += 1 + e = reinterpret(UInt16,xi) & exponent_mask + H[exponent_bias-Int(e >> 10)] += 1 end # test that the most frequent exponents occur at 50%, 25%, 12.5% etc. From ca7aef988f4e43c12dbca8487c024a14dd28a93f Mon Sep 17 00:00:00 2001 From: Milan Date: Fri, 19 Feb 2021 18:34:01 +0000 Subject: [PATCH 3/3] type constraints to avoid segfault --- src/randfloat.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/randfloat.jl b/src/randfloat.jl index 7830e7e..383d9c4 100644 --- a/src/randfloat.jl +++ b/src/randfloat.jl @@ -84,7 +84,7 @@ function randfloat(rng::Random.AbstractRNG,::Type{Float16}) end # use stdlib default RNG as a default here too -randfloat(::Type{T}=Float64) where T = randfloat(GLOBAL_RNG,T) +randfloat(::Type{T}=Float64) where {T<:Base.IEEEFloat} = randfloat(GLOBAL_RNG,T) randfloat(rng::Random.AbstractRNG) = randfloat(rng,Float64) # randfloat for arrays - in-place @@ -96,7 +96,7 @@ function randfloat!(rng::Random.AbstractRNG, A::AbstractArray{T}) where T end # randfloat for arrays with memory allocation -randfloat(rng::Random.AbstractRNG, ::Type{T}, dims::Integer...) where T = randfloat!(rng, Array{T}(undef,dims)) -randfloat(rng::Random.AbstractRNG, dims::Integer...) = randfloat!(rng, Array{Float64}(undef,dims)) -randfloat( ::Type{T}, dims::Integer...) where T = randfloat!(GLOBAL_RNG, Array{T}(undef,dims)) -randfloat( dims::Integer...) = randfloat!(GLOBAL_RNG, Array{Float64}(undef,dims)) \ No newline at end of file +randfloat(rng::Random.AbstractRNG, ::Type{T}, dims::Integer...) where {T<:Base.IEEEFloat} = randfloat!(rng, Array{T}(undef,dims)) +randfloat(rng::Random.AbstractRNG, dims::Integer...) = randfloat!(rng, Array{Float64}(undef,dims)) +randfloat(::Type{T}, dims::Integer...) where {T<:Base.IEEEFloat} = randfloat!(GLOBAL_RNG, Array{T}(undef,dims)) +randfloat( dims::Integer...) = randfloat!(GLOBAL_RNG, Array{Float64}(undef,dims)) \ No newline at end of file