From 0f2d3eac88164399f91a9ca38151af166f735c1f Mon Sep 17 00:00:00 2001 From: Rafael Fourquet Date: Fri, 24 Oct 2014 22:09:49 +0530 Subject: [PATCH 1/2] dSFMT: use fill_array_* API instead of genrand_* API Both APIs are not very compatible, so let's use fill_array_* functions, which are faster but require MersenneTwister to have an additional array field. --- base/dSFMT.jl | 19 ++++++++++++++ base/random.jl | 69 +++++++++++++++++++++++++++++++++++++++++++------- 2 files changed, 79 insertions(+), 9 deletions(-) diff --git a/base/dSFMT.jl b/base/dSFMT.jl index bddf6e2ee06b5..a51e308f54ae8 100644 --- a/base/dSFMT.jl +++ b/base/dSFMT.jl @@ -6,6 +6,7 @@ export DSFMT_state, dsfmt_get_min_array_size, dsfmt_get_idstring, dsfmt_genrand_close1_open2, dsfmt_gv_genrand_close1_open2, dsfmt_genrand_close_open, dsfmt_gv_genrand_close_open, dsfmt_genrand_uint32, dsfmt_gv_genrand_uint32, + dsfmt_fill_array_close_open!, dsfmt_fill_array_close1_open2!, win32_SystemFunction036! type DSFMT_state @@ -95,6 +96,24 @@ function dsfmt_gv_genrand_uint32() ()) end +# precondition for dsfmt_fill_array_*: +# the underlying C array must be 16-byte aligned, which is the case for "Array" +function dsfmt_fill_array_close1_open2!(s::DSFMT_state, A::Array{Float64}, n::Int) + @assert dsfmt_min_array_size <= n <= length(A) && iseven(n) + ccall((:dsfmt_fill_array_close1_open2,:libdSFMT), + Void, + (Ptr{Void}, Ptr{Float64}, Int), + s.val, A, n) +end + +function dsfmt_fill_array_close_open!(s::DSFMT_state, A::Array{Float64}, n::Int) + @assert dsfmt_min_array_size <= n <= length(A) && iseven(n) + ccall((:dsfmt_fill_array_close_open,:libdSFMT), + Void, + (Ptr{Void}, Ptr{Float64}, Int), + s.val, A, n) +end + ## Windows entropy @windows_only begin diff --git a/base/random.jl b/base/random.jl index 909c031ba2787..c40262125bc3e 100644 --- a/base/random.jl +++ b/base/random.jl @@ -14,19 +14,39 @@ abstract AbstractRNG type MersenneTwister <: AbstractRNG state::DSFMT_state seed::Union(Uint32,Vector{Uint32}) + vals::Vector{Float64} + idx::Int function MersenneTwister(seed::Vector{Uint32}) state = DSFMT_state() dsfmt_init_by_array(state, seed) - return new(state, seed) + return new(state, seed, Array(Float64, dsfmt_get_min_array_size()), dsfmt_get_min_array_size()) end MersenneTwister(seed=0) = MersenneTwister(make_seed(seed)) end +## Low level API for MersenneTwister + +function gen_rand(r::MersenneTwister) + dsfmt_fill_array_close1_open2!(r.state, r.vals, length(r.vals)) + r.idx = 0 +end + +@inline gen_rand_maybe(r::MersenneTwister) = r.idx == length(r.vals) && 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]) +@inline rand_inbounds(r::MersenneTwister) = rand_close1_open2_inbounds(r) - 1.0 + +# produce Float64 values +@inline rand_close1_open2(r::MersenneTwister) = (gen_rand_maybe(r); rand_close1_open2_inbounds(r)) +@inline rand(r::MersenneTwister) = (gen_rand_maybe(r); rand_inbounds(r)) + function srand(r::MersenneTwister, seed) r.seed = seed dsfmt_init_gen_rand(r.state, seed) + r.idx = length(r.vals) return r end @@ -96,9 +116,6 @@ rand(::Type{Float16}) = float16(rand()) rand{T<:Real}(::Type{Complex{T}}) = complex(rand(T),rand(T)) - -rand(r::MersenneTwister) = dsfmt_genrand_close_open(r.state) - ## random integers dsfmt_randui32() = dsfmt_gv_genrand_uint32() @@ -142,6 +159,41 @@ function rand!{T}(r::AbstractRNG, A::AbstractArray{T}) A end +function rand_AbstractArray_Float64!(r::MersenneTwister, A::AbstractArray{Float64}) + n = length(A) + # what follows is equivalent to this simple loop but more efficient: + # for i=1:n + # @inbounds A[i] = rand(r) + # end + m = 0 + while m < n + s = length(r.vals) - r.idx + if s == 0 + gen_rand(r) + s = length(r.vals) + end + m2 = min(n, m+s) + for i=m+1:m2 + @inbounds A[i] = rand_inbounds(r) + end + m = m2 + end + A +end + +rand!(r::MersenneTwister, A::AbstractArray{Float64}) = rand_AbstractArray_Float64!(r, A) + +function rand!(r::MersenneTwister, A::Array{Float64}) + n = length(A) + if n < dsfmt_get_min_array_size() + rand_AbstractArray_Float64!(r, A) + else + dsfmt_fill_array_close_open!(r.state, A, 2*(n รท 2)) + isodd(n) && (A[n] = rand(r)) + end + A +end + rand(T::Type, dims::Dims) = rand!(Array(T, dims)) rand{T<:Number}(::Type{T}) = error("no random number generator for type $T; try a more specific type") rand{T<:Number}(::Type{T}, dims::Int...) = rand(T, dims) @@ -737,11 +789,10 @@ ziggurat_nor_r = 3.6541528853610087963519472518 ziggurat_nor_inv_r = inv(ziggurat_nor_r) ziggurat_exp_r = 7.6971174701310497140446280481 -rand(state::DSFMT_state) = dsfmt_genrand_close_open(state) randi() = reinterpret(Uint64,dsfmt_gv_genrand_close1_open2()) & 0x000fffffffffffff -randi(state::DSFMT_state) = reinterpret(Uint64,dsfmt_genrand_close1_open2(state)) & 0x000fffffffffffff +@inline randi(rng::MersenneTwister) = reinterpret(Uint64, rand_close1_open2(rng)) & 0x000fffffffffffff for (lhs, rhs) in (([], []), - ([:(state::DSFMT_state)], [:state])) + ([:(rng::MersenneTwister)], [:rng])) @eval begin function randmtzig_randn($(lhs...)) @inbounds begin @@ -787,9 +838,9 @@ for (lhs, rhs) in (([], []), end randn() = randmtzig_randn() -randn(rng::MersenneTwister) = randmtzig_randn(rng.state) +randn(rng::MersenneTwister) = randmtzig_randn(rng) randn!(A::Array{Float64}) = (for i = 1:length(A);A[i] = randmtzig_randn();end;A) -randn!(rng::MersenneTwister, A::Array{Float64}) = (for i = 1:length(A);A[i] = randmtzig_randn(rng.state);end;A) +randn!(rng::MersenneTwister, A::Array{Float64}) = (for i = 1:length(A);A[i] = randmtzig_randn(rng);end;A) randn(dims::Dims) = randn!(Array(Float64, dims)) randn(dims::Int...) = randn!(Array(Float64, dims...)) randn(rng::MersenneTwister, dims::Dims) = randn!(rng, Array(Float64, dims)) From 3f6bd03c74bff8ddd809990bf55603489d087b3f Mon Sep 17 00:00:00 2001 From: Rafael Fourquet Date: Mon, 27 Oct 2014 20:04:41 +0530 Subject: [PATCH 2/2] make Julia global RNG be an instance of MersenneTwister The state of the global RNG was previously handled by libdSFMT. This commit allows the global RNG to benefit speed improvements resulting from the use of fill_array_* functions (cf. previous commit). However, scalar calls to rand() (i.e. producing one value instead of filling an array) using the global RNG are now slower, as well as those filling non-Float64 arrays (less than twice as slow though). --- base/random.jl | 37 ++++++++++++++++++++++++------------- 1 file changed, 24 insertions(+), 13 deletions(-) diff --git a/base/random.jl b/base/random.jl index c40262125bc3e..b1fb899c964d5 100644 --- a/base/random.jl +++ b/base/random.jl @@ -41,7 +41,11 @@ end # produce Float64 values @inline rand_close1_open2(r::MersenneTwister) = (gen_rand_maybe(r); rand_close1_open2_inbounds(r)) -@inline rand(r::MersenneTwister) = (gen_rand_maybe(r); rand_inbounds(r)) +@inline rand_close_open(r::MersenneTwister) = (gen_rand_maybe(r); rand_inbounds(r)) + +# this is similar to `dsfmt_genrand_uint32` from dSFMT.h: +@inline rand_ui32(r::MersenneTwister) = reinterpret(Uint64, rand_close1_open2(r)) % Uint32 + function srand(r::MersenneTwister, seed) r.seed = seed @@ -80,8 +84,10 @@ __init__() = srand() ## srand() function srand(seed::Vector{Uint32}) - global RANDOM_SEED = seed - dsfmt_gv_init_by_array(seed) + GLOBAL_RNG.seed = seed + dsfmt_init_by_array(GLOBAL_RNG.state, seed) + GLOBAL_RNG.idx = length(GLOBAL_RNG.vals) + return GLOBAL_RNG end srand(n::Integer) = srand(make_seed(n)) @@ -106,10 +112,16 @@ function srand(filename::String, n::Integer) end srand(filename::String) = srand(filename, 4) +## Global RNG + +const GLOBAL_RNG = MersenneTwister() +globalRNG() = GLOBAL_RNG + ## random floating point values -rand(::Type{Float64}) = dsfmt_gv_genrand_close_open() -rand() = dsfmt_gv_genrand_close_open() +rand(r::MersenneTwister=GLOBAL_RNG) = rand_close_open(r) + +rand(::Type{Float64}) = rand() rand(::Type{Float32}) = float32(rand()) rand(::Type{Float16}) = float16(rand()) @@ -118,13 +130,10 @@ rand{T<:Real}(::Type{Complex{T}}) = complex(rand(T),rand(T)) ## random integers -dsfmt_randui32() = dsfmt_gv_genrand_uint32() -dsfmt_randui64() = uint64(dsfmt_randui32()) | (uint64(dsfmt_randui32())<<32) - rand(::Type{Uint8}) = rand(Uint32) % Uint8 rand(::Type{Uint16}) = rand(Uint32) % Uint16 -rand(::Type{Uint32}) = dsfmt_randui32() -rand(::Type{Uint64}) = dsfmt_randui64() +rand(::Type{Uint32}) = rand_ui32(GLOBAL_RNG) +rand(::Type{Uint64}) = uint64(rand(Uint32)) <<32 | rand(Uint32) rand(::Type{Uint128}) = uint128(rand(Uint64))<<64 | rand(Uint64) rand(::Type{Int8}) = rand(Uint32) % Int8 @@ -194,6 +203,9 @@ function rand!(r::MersenneTwister, A::Array{Float64}) A end +rand!(A::AbstractArray{Float64}) = rand!(GLOBAL_RNG, A) +rand!(A::Array{Float64}) = rand!(GLOBAL_RNG, A) + rand(T::Type, dims::Dims) = rand!(Array(T, dims)) rand{T<:Number}(::Type{T}) = error("no random number generator for type $T; try a more specific type") rand{T<:Number}(::Type{T}, dims::Int...) = rand(T, dims) @@ -293,7 +305,7 @@ rand!(B::BitArray) = Base.bitarray_rand_fill!(B) randbool(dims::Dims) = rand!(BitArray(dims)) randbool(dims::Int...) = rand!(BitArray(dims)) -randbool() = ((dsfmt_randui32() & 1) == 1) +randbool() = ((rand(Uint32) & 1) == 1) rand(::Type{Bool}) = randbool() ## randn() - Normally distributed random numbers using Ziggurat algorithm @@ -789,8 +801,7 @@ ziggurat_nor_r = 3.6541528853610087963519472518 ziggurat_nor_inv_r = inv(ziggurat_nor_r) ziggurat_exp_r = 7.6971174701310497140446280481 -randi() = reinterpret(Uint64,dsfmt_gv_genrand_close1_open2()) & 0x000fffffffffffff -@inline randi(rng::MersenneTwister) = reinterpret(Uint64, rand_close1_open2(rng)) & 0x000fffffffffffff +@inline randi(rng::MersenneTwister=GLOBAL_RNG) = reinterpret(Uint64, rand_close1_open2(rng)) & 0x000fffffffffffff for (lhs, rhs) in (([], []), ([:(rng::MersenneTwister)], [:rng])) @eval begin