Skip to content

randfloat(Float16) #75

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
sunoru merged 3 commits into from
Feb 23, 2021
Merged

randfloat(Float16) #75

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
3a036ef
randfloat(Float16)
milankl Feb 18, 2021
52686cc
Base.exponent_mask/bias fails on 1.3
milankl Feb 18, 2021
ca7aef9
type constraints to avoid segfault
milankl Feb 19, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 29 additions & 5 deletions src/randfloat.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,8 +59,32 @@ 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(::Type{T}=Float64) where {T<:Base.IEEEFloat} = randfloat(GLOBAL_RNG,T)
randfloat(rng::Random.AbstractRNG) = randfloat(rng,Float64)

# randfloat for arrays - in-place
Expand All @@ -72,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))
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))
47 changes: 41 additions & 6 deletions test/randfloat.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -18,39 +21,44 @@ 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_mask = 0x7ff0000000000000
exponent_bias = 1023

# histogram for exponents
Expand All @@ -64,13 +72,16 @@ end
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_mask = 0x7f800000
exponent_bias = 127

# histogram for exponents
Expand All @@ -84,4 +95,28 @@ end
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

exponent_mask = 0x7c00
exponent_bias = 15

# histogram for exponents
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) & exponent_mask
H[exponent_bias-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