From 3ecd19d28fdbea7cacfd315cd49b9950ac617821 Mon Sep 17 00:00:00 2001 From: Rafael Fourquet Date: Sun, 9 Nov 2014 18:14:31 +0530 Subject: [PATCH 1/5] refactor MersenneTwister's low level API --- base/random.jl | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/base/random.jl b/base/random.jl index 599db379fbde4..b16bf5f556701 100644 --- a/base/random.jl +++ b/base/random.jl @@ -11,28 +11,36 @@ export srand, abstract AbstractRNG +const MTCacheLength = dsfmt_get_min_array_size() + type MersenneTwister <: AbstractRNG state::DSFMT_state vals::Vector{Float64} idx::Int seed::Vector{UInt32} - MersenneTwister(seed) = srand(new(DSFMT_state(), Array(Float64, dsfmt_get_min_array_size())), + MersenneTwister(seed) = srand(new(DSFMT_state(), Array(Float64, MTCacheLength)), seed) MersenneTwister() = MersenneTwister(0) end ## Low level API for MersenneTwister +@inline mt_avail(r::MersenneTwister) = MTCacheLength - r.idx +@inline mt_empty(r::MersenneTwister) = r.idx == MTCacheLength +@inline mt_setfull!(r::MersenneTwister) = r.idx = 0 +@inline mt_setempty!(r::MersenneTwister) = r.idx = MTCacheLength +@inline mt_pop!(r::MersenneTwister) = @inbounds return r.vals[r.idx+=1] + function gen_rand(r::MersenneTwister) dsfmt_fill_array_close1_open2!(r.state, r.vals, length(r.vals)) - r.idx = 0 + mt_setfull!(r) end -@inline gen_rand_maybe(r::MersenneTwister) = r.idx == length(r.vals) && gen_rand(r) +@inline gen_rand_maybe(r::MersenneTwister) = mt_empty(r) && gen_rand(r) -# precondition: r.idx < length(r.vals) -@inline rand_close1_open2_inbounds(r::MersenneTwister) = (r.idx += 1; @inbounds return r.vals[r.idx]) +# precondition: !mt_empty(r) +@inline rand_close1_open2_inbounds(r::MersenneTwister) = mt_pop!(r) @inline rand_inbounds(r::MersenneTwister) = rand_close1_open2_inbounds(r) - 1.0 # produce Float64 values @@ -46,7 +54,7 @@ end function srand(r::MersenneTwister, seed::Vector{UInt32}) r.seed = seed dsfmt_init_by_array(r.state, r.seed) - r.idx = length(r.vals) + mt_setempty!(r) return r end @@ -189,10 +197,10 @@ function rand_AbstractArray_Float64!(r::MersenneTwister, A::AbstractArray{Float6 # end m = 0 while m < n - s = length(r.vals) - r.idx + s = mt_avail(r) if s == 0 gen_rand(r) - s = length(r.vals) + s = mt_avail(r) end m2 = min(n, m+s) for i=m+1:m2 From d04f8553123a424cbcdccc93e2f3d833dcd38ee2 Mon Sep 17 00:00:00 2001 From: Rafael Fourquet Date: Sun, 9 Nov 2014 20:33:19 +0530 Subject: [PATCH 2/5] make the Float64 interval a type parameter of rand --- base/random.jl | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/base/random.jl b/base/random.jl index b16bf5f556701..4abcdd2e98932 100644 --- a/base/random.jl +++ b/base/random.jl @@ -39,16 +39,20 @@ end @inline gen_rand_maybe(r::MersenneTwister) = mt_empty(r) && gen_rand(r) +abstract FloatInterval +type CloseOpen <: FloatInterval end +type Close1Open2 <: FloatInterval end + # precondition: !mt_empty(r) -@inline rand_close1_open2_inbounds(r::MersenneTwister) = mt_pop!(r) -@inline rand_inbounds(r::MersenneTwister) = rand_close1_open2_inbounds(r) - 1.0 +@inline rand_inbounds(r::MersenneTwister, ::Type{Close1Open2}) = mt_pop!(r) +@inline rand_inbounds(r::MersenneTwister, ::Type{CloseOpen}) = rand_inbounds(r, Close1Open2) - 1.0 +@inline rand_inbounds(r::MersenneTwister) = rand_inbounds(r, CloseOpen) # produce Float64 values -@inline rand_close1_open2(r::MersenneTwister) = (gen_rand_maybe(r); rand_close1_open2_inbounds(r)) -@inline rand_close_open(r::MersenneTwister) = (gen_rand_maybe(r); rand_inbounds(r)) +@inline rand{I<:FloatInterval}(r::MersenneTwister, ::Type{I}) = (gen_rand_maybe(r); rand_inbounds(r, I)) # this is similar to `dsfmt_genrand_uint32` from dSFMT.h: -@inline rand_ui32(r::MersenneTwister) = reinterpret(UInt64, rand_close1_open2(r)) % UInt32 +@inline rand_ui32(r::MersenneTwister) = reinterpret(UInt64, rand(r, Close1Open2)) % UInt32 function srand(r::MersenneTwister, seed::Vector{UInt32}) @@ -134,7 +138,7 @@ globalRNG() = GLOBAL_RNG # rand: a non-specified RNG defaults to GLOBAL_RNG -@inline rand() = rand_close_open(GLOBAL_RNG) +@inline rand() = rand(GLOBAL_RNG, CloseOpen) @inline rand(T::Type) = rand(GLOBAL_RNG, T) rand(::()) = rand(GLOBAL_RNG, ()) # needed to resolve ambiguity rand(dims::Dims) = rand(GLOBAL_RNG, dims) @@ -145,10 +149,10 @@ rand!(A::AbstractArray) = rand!(GLOBAL_RNG, A) ## random floating point values -@inline rand(r::AbstractRNG) = rand_close_open(r) +@inline rand(r::AbstractRNG) = rand(r, CloseOpen) # MersenneTwister -rand(r::MersenneTwister, ::Type{Float64}) = rand_close_open(r) +rand(r::MersenneTwister, ::Type{Float64}) = rand(r, CloseOpen) rand{T<:Union(Float16, Float32)}(r::MersenneTwister, ::Type{T}) = convert(T, rand(r, Float64)) ## random integers (MersenneTwister) @@ -819,7 +823,7 @@ const ziggurat_nor_r = 3.6541528853610087963519472518 const ziggurat_nor_inv_r = inv(ziggurat_nor_r) const ziggurat_exp_r = 7.6971174701310497140446280481 -@inline randi(rng::MersenneTwister=GLOBAL_RNG) = reinterpret(Uint64, rand_close1_open2(rng)) & 0x000fffffffffffff +@inline randi(rng::MersenneTwister=GLOBAL_RNG) = reinterpret(Uint64, rand(rng, Close1Open2)) & 0x000fffffffffffff function randmtzig_randn(rng::MersenneTwister=GLOBAL_RNG) @inbounds begin From d3eb26bd87e935081124472ac875d4280a4921e3 Mon Sep 17 00:00:00 2001 From: Rafael Fourquet Date: Sun, 9 Nov 2014 22:42:02 +0530 Subject: [PATCH 3/5] implement faster rand!{T<:IntTypes}(::MersenneTwister, ::Array{T}) --- base/random.jl | 56 ++++++++++++++++++++++++++++++++++++++++++-------- test/random.jl | 4 ++-- 2 files changed, 49 insertions(+), 11 deletions(-) diff --git a/base/random.jl b/base/random.jl index 4abcdd2e98932..6ba211327671c 100644 --- a/base/random.jl +++ b/base/random.jl @@ -54,6 +54,8 @@ type Close1Open2 <: FloatInterval end # this is similar to `dsfmt_genrand_uint32` from dSFMT.h: @inline rand_ui32(r::MersenneTwister) = reinterpret(UInt64, rand(r, Close1Open2)) % UInt32 +@inline rand_ui52_raw(r::MersenneTwister) = reinterpret(UInt64, rand(r, Close1Open2)) +@inline rand_ui2x52_raw(r::MersenneTwister) = (((rand_ui52_raw(r) % UInt128) << 64) | rand_ui52_raw(r)) function srand(r::MersenneTwister, seed::Vector{UInt32}) r.seed = seed @@ -193,11 +195,10 @@ end # MersenneTwister -function rand_AbstractArray_Float64!(r::MersenneTwister, A::AbstractArray{Float64}) - n = length(A) +function rand_AbstractArray_Float64!{I<:FloatInterval}(r::MersenneTwister, A::AbstractArray{Float64}, n=length(A), ::Type{I}=CloseOpen) # what follows is equivalent to this simple loop but more efficient: # for i=1:n - # @inbounds A[i] = rand(r) + # @inbounds A[i] = rand(r, I) # end m = 0 while m < n @@ -208,7 +209,7 @@ function rand_AbstractArray_Float64!(r::MersenneTwister, A::AbstractArray{Float6 end m2 = min(n, m+s) for i=m+1:m2 - @inbounds A[i] = rand_inbounds(r) + @inbounds A[i] = rand_inbounds(r, I) end m = m2 end @@ -217,13 +218,50 @@ end rand!(r::MersenneTwister, A::AbstractArray{Float64}) = rand_AbstractArray_Float64!(r, A) -function rand!(r::MersenneTwister, A::Array{Float64}) - n = length(A) +fill_array!(s::DSFMT_state, A::Array{Float64}, n::Int, ::Type{CloseOpen}) = dsfmt_fill_array_close_open!(s, A, n) +fill_array!(s::DSFMT_state, A::Array{Float64}, n::Int, ::Type{Close1Open2}) = dsfmt_fill_array_close1_open2!(s, A, n) + +function rand!{I<:FloatInterval}(r::MersenneTwister, A::Array{Float64}, n=length(A), ::Type{I}=CloseOpen) if n < dsfmt_get_min_array_size() - rand_AbstractArray_Float64!(r, A) + rand_AbstractArray_Float64!(r, A, n, I) else - dsfmt_fill_array_close_open!(r.state, A, 2*(n ÷ 2)) - isodd(n) && (A[n] = rand(r)) + fill_array!(r.state, A, 2*(n ÷ 2), I) + isodd(n) && (A[n] = rand(r, I)) + end + A +end + +function rand!(r::MersenneTwister, A::Array{UInt128}, n=length(A)) + Af = pointer_to_array(convert(Ptr{Float64}, pointer(A)), 2n) + i = n + while true + rand!(r, Af, 2i, Close1Open2) + n < 5 && break + i = 0 + @inbounds while n-i >= 5 + u = A[i+=1] + A[n] $= u << 48 + A[n-=1] $= u << 36 + A[n-=1] $= u << 24 + A[n-=1] $= u << 12 + n-=1 + end + end + if n > 0 + u = rand_ui2x52_raw(r) + for i = 1:n + @inbounds A[i] $= u << 12*i + end + end + A +end + +function rand!{T<:Union(Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Int128)}(r::MersenneTwister, A::Array{T}) + n=length(A) + n128 = n * sizeof(T) ÷ 16 + rand!(r, pointer_to_array(convert(Ptr{UInt128}, pointer(A)), n128)) + for i = 16*n128÷sizeof(T)+1:n + @inbounds A[i] = rand(r, T) end A end diff --git a/test/random.jl b/test/random.jl index 9271c49981087..a53cbd67a2ba6 100644 --- a/test/random.jl +++ b/test/random.jl @@ -28,8 +28,8 @@ rand!(MersenneTwister(0), A) @test rand(MersenneTwister(0), Int64, 1) == [172014471070449746] A = zeros(Int64, 2, 2) rand!(MersenneTwister(0), A) -@test A == [ 172014471070449746 -193283627354378518; - -4679130500738884555 -9008350441255501549] +@test A == [858542123778948672 5715075217119798169; + 8690327730555225005 8435109092665372532] # randn @test randn(MersenneTwister(42)) == -0.5560268761463861 From f6f0a8269de9e3fb39832c2d78056d713d130e60 Mon Sep 17 00:00:00 2001 From: Rafael Fourquet Date: Mon, 10 Nov 2014 03:31:26 +0530 Subject: [PATCH 4/5] implement faster rand!{T<:Union(Float16,Float32)}(::MersenneTwister, ::Array{T}) --- base/random.jl | 39 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/base/random.jl b/base/random.jl index 6ba211327671c..07f1b833d7d6e 100644 --- a/base/random.jl +++ b/base/random.jl @@ -231,6 +231,45 @@ function rand!{I<:FloatInterval}(r::MersenneTwister, A::Array{Float64}, n=length A end +@inline mask128(u::UInt128, ::Type{Float16}) = (u & 0x03ff03ff03ff03ff03ff03ff03ff03ff) | 0x3c003c003c003c003c003c003c003c00 +@inline mask128(u::UInt128, ::Type{Float32}) = (u & 0x007fffff007fffff007fffff007fffff) | 0x3f8000003f8000003f8000003f800000 + +function rand!{T<:Union(Float16, Float32)}(r::MersenneTwister, A::Array{T}, ::Type{Close1Open2}) + n = length(A) + n128 = n * sizeof(T) ÷ 16 + rand!(r, pointer_to_array(convert(Ptr{Float64}, pointer(A)), 2*n128), 2*n128, Close1Open2) + A128 = pointer_to_array(convert(Ptr{UInt128}, pointer(A)), n128) + @inbounds for i in 1:n128 + u = A128[i] + u $= u << 26 + # at this point, the 64 low bits of u, "k" being the k-th bit of A128[i] and "+" the bit xor, are: + # [..., 58+32,..., 53+27, 52+26, ..., 33+7, 32+6, ..., 27+1, 26, ..., 1] + # the bits needing to be random are + # [1:10, 17:26, 33:42, 49:58] (for Float16) + # [1:23, 33:55] (for Float32) + # this is obviously satisfied on the 32 low bits side, and on the high side, the entropy comes + # from bits 33:52 of A128[i] and then from bits 27:32 (which are discarded on the low side) + # this is similar for the 64 high bits of u + A128[i] = mask128(u, T) + end + for i in 16*n128÷sizeof(T)+1:n + @inbounds A[i] = rand(r, T) + one(T) + end + A +end + +function rand!{T<:Union(Float16, Float32)}(r::MersenneTwister, A::Array{T}, ::Type{CloseOpen}) + rand!(r, A, Close1Open2) + I32 = one(Float32) + for i in 1:length(A) + @inbounds A[i] = T(Float32(A[i])-I32) # faster than "A[i] -= one(T)" for T==Float16 + end + A +end + +rand!{T<:Union(Float16, Float32)}(r::MersenneTwister, A::Array{T}) = rand!(r, A, CloseOpen) + + function rand!(r::MersenneTwister, A::Array{UInt128}, n=length(A)) Af = pointer_to_array(convert(Ptr{Float64}, pointer(A)), 2n) i = n From 912e9adb0f8bdeb01da99bb54ab9f0b103be938a Mon Sep 17 00:00:00 2001 From: Rafael Fourquet Date: Tue, 11 Nov 2014 17:28:52 +0530 Subject: [PATCH 5/5] add rand! tests --- test/random.jl | 40 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/test/random.jl b/test/random.jl index a53cbd67a2ba6..ab871c2e7f4fe 100644 --- a/test/random.jl +++ b/test/random.jl @@ -182,3 +182,43 @@ i8257 = 1:1/3:100 for i = 1:100 @test rand(i8257) in i8257 end + +# test code paths of rand! + +let mt = MersenneTwister(0) + A128 = Array(UInt128, 0) + @test length(rand!(mt, A128)) == 0 + for (i,n) in enumerate([1, 3, 5, 6, 10, 11, 30]) + resize!(A128, n) + rand!(mt, A128) + @test length(A128) == n + @test A128[end] == UInt128[0x15de6b23025813ad129841f537a04e40, + 0xcfa4db38a2c65bc4f18c07dc91125edf, + 0x33bec08136f19b54290982449b3900d5, + 0xde41af3463e74cb830dad4add353ca20, + 0x066d8695ebf85f833427c93416193e1f, + 0x48fab49cc9fcee1c920d6dae629af446, + 0x4b54632b4619f4eca22675166784d229][i] + + end + + srand(mt,0) + for (i,T) in enumerate([Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Int128, Float16, Float32]) + A = Array(T, 16) + B = Array(T, 31) + rand!(mt, A) + rand!(mt, B) + @test A[end] == Any[21,0x7b,17385,0x3086,-1574090021,0xadcb4460,6797283068698303107,0x4e91c9c4d4f5f759, + -3482609696641744459568613291754091152,float16(0.03125),0.68733835f0][i] + + @test B[end] == Any[49,0x65,-3725,0x719d,814246081,0xdf61843a,-1603010949539670188,0x5e4ca1658810985d, + -33032345278809823492812856023466859769,float16(0.9346),0.5929704f0][i] + end + + srand(mt,0) + AF64 = Array(Float64, Base.Random.dsfmt_get_min_array_size()-1) + @test rand!(mt, AF64)[end] == 0.957735065345398 + @test rand!(mt, AF64)[end] == 0.6492481059865669 + resize!(AF64, 2*length(mt.vals)) + @test Base.Random.rand_AbstractArray_Float64!(mt, AF64)[end] == 0.432757268470779 +end