From b8bc81677ea4075d161a464dbb33f0eef33aa291 Mon Sep 17 00:00:00 2001 From: Rafael Fourquet Date: Sat, 17 Oct 2020 18:45:43 +0200 Subject: [PATCH] quite faster rand(::MersenneTwister, ::Type{Int}) ... (#37914) ... when the quest for superficial beauty leads to performance ... The recent implementation of `show` for `MersenneTwister` was not ideal, as a number of book-keeping variables had to be introduced as fields of MT; in particular, as generating the cache for ints (using the generic `rand!` for integer arrays) was consuming random `Float64` numbers, 4 integers had to be shown in `show` only to reproduce the state of the ints cache. In total 8 integers were shown. But it is not so difficult to improve a bit, thanks to two features of the internal cache: 1) it's 16-byte aligned, so the dSFMT low-level routine can be called directly on it (whereas the generic `rand!` for integers has to take care that the same stream is produced whatever the alignment) 2) it can be resized: dSFMT randomizing only 52 out of 64 bits, i.e. a bit more than 80% of the bits, the trick is to count the total number of needed bits, grow the array to a size such that dSFMT produces these needed bits, and then condense these bits back into a 100% randomized array of the original size As a net result, two variables could be deleted which `show` doesn't needs to display anymore. Purely as a side effect, scalar generation of `Int64`/`UInt64` has a speedup of about 1.8x, and about 1.6x for `Int128`/`UInt128`, and about 1.3x or 1.5x for generation in a range of integers of size with less than 64 bits, e.g. `rand(1:9)` (at least on this machine...) --- stdlib/Random/docs/src/index.md | 18 ++--- stdlib/Random/src/RNGs.jl | 126 +++++++++++++++++--------------- stdlib/Random/src/misc.jl | 14 ++-- stdlib/Random/test/runtests.jl | 36 +++++---- stdlib/UUIDs/src/UUIDs.jl | 6 +- 5 files changed, 108 insertions(+), 92 deletions(-) diff --git a/stdlib/Random/docs/src/index.md b/stdlib/Random/docs/src/index.md index 0de918f10a1af..11214ca733ff6 100644 --- a/stdlib/Random/docs/src/index.md +++ b/stdlib/Random/docs/src/index.md @@ -145,22 +145,22 @@ Scalar and array methods for `Die` now work as expected: ```jldoctest Die; setup = :(Random.seed!(1)) julia> rand(Die) -Die(10) +Die(15) julia> rand(MersenneTwister(0), Die) -Die(16) +Die(11) julia> rand(Die, 3) 3-element Vector{Die}: + Die(18) Die(5) - Die(20) - Die(9) + Die(4) julia> a = Vector{Die}(undef, 3); rand!(a) 3-element Vector{Die}: - Die(11) + Die(5) Die(20) - Die(10) + Die(15) ``` #### A simple sampler without pre-computed data @@ -173,13 +173,13 @@ In order to define random generation out of objects of type `S`, the following m julia> Random.rand(rng::AbstractRNG, d::Random.SamplerTrivial{Die}) = rand(rng, 1:d[].nsides); julia> rand(Die(4)) -2 +3 julia> rand(Die(4), 3) 3-element Vector{Any}: - 1 4 - 2 + 1 + 1 ``` Given a collection type `S`, it's currently assumed that if `rand(::S)` is defined, an object of type `eltype(S)` will be produced. In the last example, a `Vector{Any}` is produced; the reason is that `eltype(Die) == Any`. The remedy is to define `Base.eltype(::Type{Die}) = Int`. diff --git a/stdlib/Random/src/RNGs.jl b/stdlib/Random/src/RNGs.jl index 6345fbfc23811..7bc0a1c242ff4 100644 --- a/stdlib/Random/src/RNGs.jl +++ b/stdlib/Random/src/RNGs.jl @@ -90,11 +90,9 @@ mutable struct MersenneTwister <: AbstractRNG adv_jump::BigInt # number of skipped Float64 values via randjump adv_vals::Int64 # state of advance when vals is filled-up adv_ints::Int64 # state of advance when ints is filled-up - adv_vals_pre::Int64 # state of advance when vals is filled-up before ints - adv_idxF_pre::Int # value of idxF before ints is filled-up function MersenneTwister(seed, state, vals, ints, idxF, idxI, - adv, adv_jump, adv_vals, adv_ints, adv_vals_pre, adv_idxF_pre) + adv, adv_jump, adv_vals, adv_ints) length(vals) == MT_CACHE_F && 0 <= idxF <= MT_CACHE_F || throw(DomainError((length(vals), idxF), "`length(vals)` and `idxF` must be consistent with $MT_CACHE_F")) @@ -102,7 +100,7 @@ mutable struct MersenneTwister <: AbstractRNG throw(DomainError((length(ints), idxI), "`length(ints)` and `idxI` must be consistent with $MT_CACHE_I")) new(seed, state, vals, ints, idxF, idxI, - adv, adv_jump, adv_vals, adv_ints, adv_vals_pre, adv_idxF_pre) + adv, adv_jump, adv_vals, adv_ints) end end @@ -110,7 +108,7 @@ MersenneTwister(seed::Vector{UInt32}, state::DSFMT_state) = MersenneTwister(seed, state, Vector{Float64}(undef, MT_CACHE_F), Vector{UInt128}(undef, MT_CACHE_I >> 4), - MT_CACHE_F, 0, 0, 0, -1, -1, -1, -1) + MT_CACHE_F, 0, 0, 0, -1, -1) """ MersenneTwister(seed) @@ -161,15 +159,12 @@ function copy!(dst::MersenneTwister, src::MersenneTwister) dst.adv_jump = src.adv_jump dst.adv_vals = src.adv_vals dst.adv_ints = src.adv_ints - dst.adv_vals_pre = src.adv_vals_pre - dst.adv_idxF_pre = src.adv_idxF_pre dst end copy(src::MersenneTwister) = MersenneTwister(copy(src.seed), copy(src.state), copy(src.vals), copy(src.ints), - src.idxF, src.idxI, src.adv, src.adv_jump, src.adv_vals, src.adv_ints, - src.adv_vals_pre, src.adv_idxF_pre) + src.idxF, src.idxI, src.adv, src.adv_jump, src.adv_vals, src.adv_ints) ==(r1::MersenneTwister, r2::MersenneTwister) = @@ -191,17 +186,18 @@ function show(io::IO, rng::MersenneTwister) print(io, "MersenneTwister($seed_str, (") # state adv = Integer[rng.adv_jump, rng.adv] - if rng.adv_vals != -1 - push!(adv, rng.adv_vals, rng.idxF) + if rng.adv_vals != -1 || rng.adv_ints != -1 + if rng.adv_vals == -1 + @assert rng.idxF == MT_CACHE_F + push!(adv, 0, 0) # "(0, 0)" is nicer on the eyes than (-1, 1002) + else + push!(adv, rng.adv_vals, rng.idxF) + end end - if rng.adv_ints != -1 # then rng.adv_vals is always != -1 + if rng.adv_ints != -1 idxI = (length(rng.ints)*16 - rng.idxI) / 8 # 8 represents one Int64 idxI = Int(idxI) # idxI should always be an integer when using public APIs - push!(adv, - rng.adv_ints, - rng.adv_vals_pre == -1 ? 0 : rng.adv_vals_pre, - rng.adv_vals_pre == -1 ? 0 : rng.adv_idxF_pre, - idxI) + push!(adv, rng.adv_ints, idxI) end join(io, adv, ", ") print(io, "))") @@ -217,8 +213,6 @@ function reset_caches!(r::MersenneTwister) mt_setempty!(r, UInt128) r.adv_vals = -1 r.adv_ints = -1 - r.adv_vals_pre = -1 - r.adv_idxF_pre = -1 r end @@ -261,9 +255,33 @@ mt_avail(r::MersenneTwister, ::Type{T}) where {T<:BitInteger} = function mt_setfull!(r::MersenneTwister, ::Type{<:BitInteger}) r.adv_ints = r.adv - r.adv_vals_pre = r.adv_vals - r.adv_idxF_pre = r.idxF - rand!(r, r.ints) + ints = r.ints + + @assert length(ints) == 501 + # dSFMT natively randomizes 52 out of 64 bits of each UInt64 words, + # i.e. 12 bits are missing; + # by generating 5 words == 5*52 == 260 bits, we can fully + # randomize 4 UInt64 = 256 bits; IOW, at the array level, we must + # randomize ceil(501*1.25) = 627 UInt128 words (with 2*52 bits each), + # which we then condense into fully randomized 501 UInt128 words + + len = 501 + 126 # 126 == ceil(501 / 4) + resize!(ints, len) + p = pointer(ints) # must be *after* resize! + GC.@preserve r fill_array!(r, Ptr{Float64}(p), len*2, CloseOpen12_64()) + + k = 501 + n = 0 + @inbounds while n != 500 + u = ints[k+=1] + ints[n+=1] ⊻= u << 48 + ints[n+=1] ⊻= u << 36 + ints[n+=1] ⊻= u << 24 + ints[n+=1] ⊻= u << 12 + end + @assert k == len - 1 + @inbounds ints[501] ⊻= ints[len] << 48 + resize!(ints, 501) r.idxI = MT_CACHE_I end @@ -740,17 +758,17 @@ jump!(r::MersenneTwister, steps::Integer) = copy!(r, jump(r, steps)) # parameters in the tuples are: # 1: .adv_jump (jump steps) # 2: .adv (number of generated floats at the DSFMT_state level since seeding, besides jumps) -# 3, 4: .adv_vals, .idxF (counters to reconstruct the float chache, optional if 5-8 not shown)) -# 5-8: .adv_ints, .adv_vals_pre, .adv_idxF_pre, .idxI (counters to reconstruct the integer chache, optional) +# 3, 4: .adv_vals, .idxF (counters to reconstruct the float chache, optional if 5-6 not shown)) +# 5, 6: .adv_ints, .idxI (counters to reconstruct the integer chache, optional) -Random.MersenneTwister(seed::Union{Integer,Vector{UInt32}}, advance::NTuple{8,Integer}) = +Random.MersenneTwister(seed::Union{Integer,Vector{UInt32}}, advance::NTuple{6,Integer}) = advance!(MersenneTwister(seed), advance...) Random.MersenneTwister(seed::Union{Integer,Vector{UInt32}}, advance::NTuple{4,Integer}) = - MersenneTwister(seed, (advance..., -1, -1, -1, -1)) + MersenneTwister(seed, (advance..., 0, 0)) Random.MersenneTwister(seed::Union{Integer,Vector{UInt32}}, advance::NTuple{2,Integer}) = - MersenneTwister(seed, (advance..., 0, 0, -1, -1, -1, -1)) + MersenneTwister(seed, (advance..., 0, 0, 0, 0)) # advances raw state (per fill_array!) of r by n steps (Float64 values) function _advance_n!(r::MersenneTwister, n::Int64, work::Vector{Float64}) @@ -775,24 +793,10 @@ function _advance_to!(r::MersenneTwister, adv::Int64, work) end function _advance_F!(r::MersenneTwister, adv_vals, idxF, work) - if adv_vals == idxF == 0 - # this case happens only when integer cache was generated before float cache - # then (0, 0) is printed instead of (-1, MT_CACHE_F) which is somewhat confusing; - # in this case, nothing to do, the float cache mustn't be filled - if r.adv_vals == -1 && r.idxF == MT_CACHE_F - return - else - throw(DomainError(n, "can't advance $r to the specified state")) - end - end - if r.adv_vals != adv_vals - _advance_to!(r, adv_vals, work) - gen_rand(r) - @assert r.adv_vals == adv_vals - end # otherwise, advancing was done automatically while generating the integer cache - + _advance_to!(r, adv_vals, work) + gen_rand(r) + @assert r.adv_vals == adv_vals r.idxF = idxF - nothing end function _advance_I!(r::MersenneTwister, adv_ints, idxI, work) @@ -802,26 +806,34 @@ function _advance_I!(r::MersenneTwister, adv_ints, idxI, work) r.idxI = 16*length(r.ints) - 8*idxI end -function advance!(r::MersenneTwister, adv_jump, adv, adv_vals, idxF, - adv_ints, adv_vals_pre, adv_idxF_pre, idxI) +function advance!(r::MersenneTwister, adv_jump, adv, adv_vals, idxF, adv_ints, idxI) adv_jump = BigInt(adv_jump) - adv, adv_vals, adv_ints, adv_vals_pre = Int64.((adv, adv_vals, adv_ints, adv_vals_pre)) - idxF, adv_idxF_pre, idxI = Int.((idxF, adv_idxF_pre, idxI)) + adv, adv_vals, adv_ints = Int64.((adv, adv_vals, adv_ints)) + idxF, idxI = Int.((idxF, idxI)) ms = dsfmt_get_min_array_size() % Int work = sizehint!(Vector{Float64}(), 2ms) - jump!(r, adv_jump) - if adv_vals_pre != -1 - _advance_F!(r, adv_vals_pre, adv_idxF_pre, work) - _advance_I!(r, adv_ints, idxI, work) - @assert r.adv_vals_pre == adv_vals_pre || - r.adv_vals_pre == -1 && adv_vals_pre == 0 - @assert r.adv_idxF_pre == adv_idxF_pre || - r.adv_idxF_pre == 1002 && adv_idxF_pre == 0 + adv_jump != 0 && jump!(r, adv_jump) + advF = (adv_vals, idxF) != (0, 0) + advI = (adv_ints, idxI) != (0, 0) + if advI && advF + @assert adv_vals != adv_ints + if adv_vals < adv_ints + _advance_F!(r, adv_vals, idxF, work) + _advance_I!(r, adv_ints, idxI, work) + else + _advance_I!(r, adv_ints, idxI, work) + _advance_F!(r, adv_vals, idxF, work) + end + elseif advF + _advance_F!(r, adv_vals, idxF, work) + elseif advI + _advance_I!(r, adv_ints, idxI, work) + else + @assert adv == 0 end - _advance_F!(r, adv_vals, idxF, work) _advance_to!(r, adv, work) r end diff --git a/stdlib/Random/src/misc.jl b/stdlib/Random/src/misc.jl index 465ba27b39332..3c09f5b30d257 100644 --- a/stdlib/Random/src/misc.jl +++ b/stdlib/Random/src/misc.jl @@ -22,14 +22,14 @@ julia> rng = MersenneTwister(1234); julia> bitrand(rng, 10) 10-element BitVector: 0 - 1 - 1 - 1 - 1 + 0 + 0 0 1 0 0 + 0 + 1 1 ``` """ @@ -53,13 +53,13 @@ number generator, see [Random Numbers](@ref). # Examples ```jldoctest julia> Random.seed!(3); randstring() -"4zSHdXlw" +"Y7m62wOj" julia> randstring(MersenneTwister(3), 'a':'z', 6) -"bzlhqn" +"ocucay" julia> randstring("ACGT") -"AGGACATT" +"ATTTGCGT" ``` !!! note diff --git a/stdlib/Random/test/runtests.jl b/stdlib/Random/test/runtests.jl index bac8049de9ee1..501fa5928bd85 100644 --- a/stdlib/Random/test/runtests.jl +++ b/stdlib/Random/test/runtests.jl @@ -47,7 +47,7 @@ let A = zeros(2, 2) end let A = zeros(2, 2) @test_throws ArgumentError rand!(MersenneTwister(0), A, 5) - @test rand(MersenneTwister(0), Int64, 1) == [2118291759721269919] + @test rand(MersenneTwister(0), Int64, 1) == [-3433174948434291912] end let A = zeros(Int64, 2, 2) rand!(MersenneTwister(0), A) @@ -253,16 +253,20 @@ let mt = MersenneTwister(0) end Random.seed!(mt, 0) + Aend = Any[] + Bend = Any[] for (i,T) in enumerate([Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Int128, Float16, Float32]) A = Vector{T}(undef, 16) B = Vector{T}(undef, 31) rand!(mt, A) rand!(mt, B) - @test A[end] == Any[21, 0x7b, 17385, 0x3086, -1574090021, 0xadcb4460, 6797283068698303107, 0xc8e6453e139271f3, - 69855512850528774484795047199183096941, Float16(0.16895), 0.21086597f0][i] - @test B[end] == Any[49, 0x65, -3725, 0x719d, 814246081, 0xdf61843a, 2120308604158549401, 0xcb28c236e9c0f608, - 61881313582466480231846019869039259750, Float16(0.38672), 0.20027375f0][i] + push!(Aend, A[end]) + push!(Bend, B[end]) end + @test Aend == Any[21, 0x7b, 17385, 0x3086, -1574090021, 0xadcb4460, 6797283068698303107, 0x68a9f9865393cfd6, + 33687499368208574024854346399216845930, Float16(0.7744), 0.97259974f0] + @test Bend == Any[49, 0x65, -3725, 0x719d, 814246081, 0xdf61843a, -3433174948434291912, 0xd461716f27c91500, + -85900088726243933988214632401750448432, Float16(0.10645), 0.13879478f0] Random.seed!(mt, 0) AF64 = Vector{Float64}(undef, Random.dsfmt_get_min_array_size()-1) @@ -592,16 +596,16 @@ end @test_throws DomainError DSFMT.DSFMT_state(zeros(Int32, rand(0:DSFMT.JN32-1))) @test_throws DomainError MersenneTwister(zeros(UInt32, 1), DSFMT.DSFMT_state(), - zeros(Float64, 10), zeros(UInt128, MT_CACHE_I>>4), 0, 0, 0, 0, -1, -1, -1, -1) + zeros(Float64, 10), zeros(UInt128, MT_CACHE_I>>4), 0, 0, 0, 0, -1, -1) @test_throws DomainError MersenneTwister(zeros(UInt32, 1), DSFMT.DSFMT_state(), - zeros(Float64, MT_CACHE_F), zeros(UInt128, MT_CACHE_I>>4), -1, 0, 0, 0, -1, -1, -1, -1) + zeros(Float64, MT_CACHE_F), zeros(UInt128, MT_CACHE_I>>4), -1, 0, 0, 0, -1, -1) @test_throws DomainError MersenneTwister(zeros(UInt32, 1), DSFMT.DSFMT_state(), - zeros(Float64, MT_CACHE_F), zeros(UInt128, MT_CACHE_I>>3), 0, 0, 0, 0, -1, -1, -1, -1) + zeros(Float64, MT_CACHE_F), zeros(UInt128, MT_CACHE_I>>3), 0, 0, 0, 0, -1, -1) @test_throws DomainError MersenneTwister(zeros(UInt32, 1), DSFMT.DSFMT_state(), - zeros(Float64, MT_CACHE_F), zeros(UInt128, MT_CACHE_I>>4), 0, -1, 0, 0, -1, -1, -1, -1) + zeros(Float64, MT_CACHE_F), zeros(UInt128, MT_CACHE_I>>4), 0, -1, 0, 0, -1, -1) # seed is private to MersenneTwister let seed = rand(UInt32, 10) @@ -842,19 +846,19 @@ end @test m == MersenneTwister(123, (200000000000000000000, 0)) rand(m) @test string(m) == "MersenneTwister(123, (200000000000000000000, 1002, 0, 1))" + @test m == MersenneTwister(123, (200000000000000000000, 1002, 0, 1)) rand(m, Int64) - @test string(m) == "MersenneTwister(123, (200000000000000000000, 2002, 0, 255, 1002, 0, 1, 1))" - @test m == MersenneTwister(123, (200000000000000000000, 2002, 0, 255, 1002, 0, 1, 1)) + @test string(m) == "MersenneTwister(123, (200000000000000000000, 2256, 0, 1, 1002, 1))" + @test m == MersenneTwister(123, (200000000000000000000, 2256, 0, 1, 1002, 1)) m = MersenneTwister(0x0ecfd77f89dcd508caa37a17ebb7556b) @test string(m) == "MersenneTwister(0xecfd77f89dcd508caa37a17ebb7556b)" rand(m, Int64) - @test string(m) == "MersenneTwister(0xecfd77f89dcd508caa37a17ebb7556b, (0, 2002, 1000, 254, 0, 0, 0, 1))" - @test m == MersenneTwister(0xecfd77f89dcd508caa37a17ebb7556b, (0, 2002, 1000, 254, 0, 0, 0, 1)) + @test string(m) == "MersenneTwister(0xecfd77f89dcd508caa37a17ebb7556b, (0, 1254, 0, 0, 0, 1))" + @test m == MersenneTwister(0xecfd77f89dcd508caa37a17ebb7556b, (0, 1254, 0, 0, 0, 1)) - # test when floats advancing is done by initializing ints, and (few) floats are then generated m = MersenneTwister(0); rand(m, Int64); rand(m) - @test string(m) == "MersenneTwister(0, (0, 2002, 1000, 255, 0, 0, 0, 1))" - @test m == MersenneTwister(0, (0, 2002, 1000, 255, 0, 0, 0, 1)) + @test string(m) == "MersenneTwister(0, (0, 2256, 1254, 1, 0, 1))" + @test m == MersenneTwister(0, (0, 2256, 1254, 1, 0, 1)) end diff --git a/stdlib/UUIDs/src/UUIDs.jl b/stdlib/UUIDs/src/UUIDs.jl index 7fb83a393be00..41d5319fec24d 100644 --- a/stdlib/UUIDs/src/UUIDs.jl +++ b/stdlib/UUIDs/src/UUIDs.jl @@ -102,7 +102,7 @@ detail that may change in the future. julia> rng = MersenneTwister(1234); julia> uuid4(rng) -UUID("196f2941-2d58-45ba-9f13-43a2532b2fa8") +UUID("7a052949-c101-4ca3-9a7e-43a2532b2fa8") ``` """ function uuid4(rng::AbstractRNG=Random.RandomDevice()) @@ -126,10 +126,10 @@ as specified by RFC 4122. julia> rng = MersenneTwister(1234); julia> u4 = uuid4(rng) -UUID("196f2941-2d58-45ba-9f13-43a2532b2fa8") +UUID("7a052949-c101-4ca3-9a7e-43a2532b2fa8") julia> u5 = uuid5(u4, "julia") -UUID("b37756f8-b0c0-54cd-a466-19b3d25683bc") +UUID("086cc5bb-2461-57d8-8068-0aed7f5b5cd1") ``` """ function uuid5(ns::UUID, name::String)