diff --git a/NEWS.md b/NEWS.md index 0052496a76f3f..e5ec77ff38f5e 100644 --- a/NEWS.md +++ b/NEWS.md @@ -169,6 +169,13 @@ Standard library changes * `randn!(::MersenneTwister, ::Array{Float64})` is faster, and as a result, for a given state of the RNG, the corresponding generated numbers have changed ([#35078]). +* A new faster algorithm ("nearly division less") is used for generating random numbers + within a range ([#29240]). + As a result, the streams of generated numbers is changed. + Also, for performance, the undocumented property that, given a seed and `a, b` of type `Int`, + `rand(a:b)` produces the same stream on 32 and 64 bits architectures, is dropped. + + #### REPL diff --git a/stdlib/LinearAlgebra/test/bunchkaufman.jl b/stdlib/LinearAlgebra/test/bunchkaufman.jl index 6422887068148..5098f818f1804 100644 --- a/stdlib/LinearAlgebra/test/bunchkaufman.jl +++ b/stdlib/LinearAlgebra/test/bunchkaufman.jl @@ -12,7 +12,7 @@ n = 10 n1 = div(n, 2) n2 = 2*n1 -Random.seed!(1234321) +Random.seed!(12343210) areal = randn(n,n)/2 aimg = randn(n,n)/2 diff --git a/stdlib/Random/docs/src/index.md b/stdlib/Random/docs/src/index.md index 5a01c4665fa11..06fbfb0ea9880 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(18) +Die(10) julia> rand(MersenneTwister(0), Die) -Die(4) +Die(16) julia> rand(Die, 3) 3-element Array{Die,1}: - Die(6) - Die(11) Die(5) + Die(20) + Die(9) julia> a = Vector{Die}(undef, 3); rand!(a) 3-element Array{Die,1}: - Die(18) - Die(6) - Die(8) + Die(11) + Die(20) + Die(10) ``` #### A simple sampler without pre-computed data @@ -173,11 +173,11 @@ 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)) -3 +2 julia> rand(Die(4), 3) 3-element Array{Any,1}: - 3 + 1 4 2 ``` diff --git a/stdlib/Random/src/RNGs.jl b/stdlib/Random/src/RNGs.jl index 110fc367deb76..e95e5cc20145c 100644 --- a/stdlib/Random/src/RNGs.jl +++ b/stdlib/Random/src/RNGs.jl @@ -602,14 +602,6 @@ for T in BitInteger_types end end -#### from a range - -for T in BitInteger_types, R=(1, Inf) # eval because of ambiguity otherwise - @eval Sampler(::Type{MersenneTwister}, r::AbstractUnitRange{$T}, ::Val{$R}) = - SamplerRangeFast(r) -end - - ### randjump # Old randjump methods are deprecated, the scalar version is in the Future module. diff --git a/stdlib/Random/src/generation.jl b/stdlib/Random/src/generation.jl index caa5e3cde3f45..329620c3ab40b 100644 --- a/stdlib/Random/src/generation.jl +++ b/stdlib/Random/src/generation.jl @@ -166,13 +166,24 @@ end ### BitInteger -# there are two implemented samplers for unit ranges, which assume that Float64 (i.e. -# 52 random bits) is the native type for the RNG: -# 1) "Fast", which is the most efficient when the underlying RNG produces rand(Float64) -# "fast enough". The tradeoff is faster creation of the sampler, but more -# consumption of entropy bits -# 2) "Default" which tries to use as few entropy bits as possible, at the cost of a -# a bigger upfront price associated with the creation of the sampler +# there are three implemented samplers for unit ranges, the two first of which +# assume that Float64 (i.e. 52 random bits) is the native type for the RNG: +# 1) "Fast" (SamplerRangeFast), which is most efficient when the underlying RNG produces +# rand(Float64) "fast enough". +# The tradeoff is faster creation of the sampler, but more consumption of entropy bits. +# 2) "Slow" (SamplerRangeInt) which tries to use as few entropy bits as possible, at the +# cost of a a bigger upfront price associated with the creation of the sampler. +# This sampler is most appropriate for slower random generators. +# 3) "Nearly Division Less" (NDL) which is generally the fastest algorithm for types of size +# up to 64 bits. This is the default for these types since Julia 1.5. +# The "Fast" algorithm can be faster than NDL when the length of the range is +# less than and close to a power of 2. + +Sampler(::Type{<:AbstractRNG}, r::AbstractUnitRange{T}, + ::Repetition) where {T<:Base.BitInteger64} = SamplerRangeNDL(r) + +Sampler(::Type{<:AbstractRNG}, r::AbstractUnitRange{T}, + ::Repetition) where {T<:Union{Int128,UInt128}} = SamplerRangeFast(r) #### helper functions @@ -224,7 +235,7 @@ function rand(rng::AbstractRNG, sp::SamplerRangeFast{UInt128,T}) where T x % T + a end -#### Default +#### "Slow" / SamplerRangeInt # remainder function according to Knuth, where rem_knuth(a, 0) = a rem_knuth(a::UInt, b::UInt) = a % (b + (b == 0)) + a * (b == 0) @@ -274,10 +285,6 @@ function SamplerRangeInt(r::AbstractUnitRange{T}, ::Type{U}) where {T,U} SamplerRangeInt{T,U}(a, bw, k, mult) # overflow ok end -Sampler(::Type{<:AbstractRNG}, r::AbstractUnitRange{T}, - ::Repetition) where {T<:BitInteger} = SamplerRangeInt(r) - - rand(rng::AbstractRNG, sp::SamplerRangeInt{T,UInt32}) where {T<:BitInteger} = (unsigned(sp.a) + rem_knuth(rand(rng, LessThan(sp.u, UInt52Raw(UInt32))), sp.k)) % T @@ -295,6 +302,41 @@ function rand(rng::AbstractRNG, sp::SamplerRangeInt{T,UInt128}) where T<:BitInte return ((sp.a % UInt128) + rem_knuth(x, sp.k)) % T end +#### Nearly Division Less + +# cf. https://arxiv.org/abs/1805.10941 (algorithm 5) + +struct SamplerRangeNDL{U<:Unsigned,T} <: Sampler{T} + a::T # first element of the range + s::U # range length or zero for full range +end + +function SamplerRangeNDL(r::AbstractUnitRange{T}) where {T} + isempty(r) && throw(ArgumentError("range must be non-empty")) + a = first(r) + U = uint_sup(T) + s = (last(r) - first(r)) % unsigned(T) % U + one(U) # overflow ok + # mod(-s, s) could be put in the Sampler object for repeated calls, but + # this would be an advantage only for very big s and number of calls + SamplerRangeNDL(a, s) +end + +function rand(rng::AbstractRNG, sp::SamplerRangeNDL{U,T}) where {U,T} + s = sp.s + x = widen(rand(rng, U)) + m = x * s + l = m % U + if l < s + t = mod(-s, s) # as s is unsigned, -s is equal to 2^L - s in the paper + while l < t + x = widen(rand(rng, U)) + m = x * s + l = m % U + end + end + (s == 0 ? x : m >> (8*sizeof(U))) % T + sp.a +end + ### BigInt diff --git a/stdlib/Random/src/misc.jl b/stdlib/Random/src/misc.jl index 91b245d2f1cb0..dffaab89de3b2 100644 --- a/stdlib/Random/src/misc.jl +++ b/stdlib/Random/src/misc.jl @@ -52,14 +52,14 @@ number generator, see [Random Numbers](@ref). # Examples ```jldoctest -julia> Random.seed!(0); randstring() -"0IPrGg0J" +julia> Random.seed!(3); randstring() +"4zSHdXlw" -julia> randstring(MersenneTwister(0), 'a':'z', 6) -"aszvqk" +julia> randstring(MersenneTwister(3), 'a':'z', 6) +"bzlhqn" julia> randstring("ACGT") -"TATCGGTC" +"AGGACATT" ``` !!! note diff --git a/stdlib/Random/test/runtests.jl b/stdlib/Random/test/runtests.jl index 20f81f73e344a..b78e3ae4b8a1f 100644 --- a/stdlib/Random/test/runtests.jl +++ b/stdlib/Random/test/runtests.jl @@ -10,7 +10,7 @@ using .Main.OffsetArrays using Random using Random.DSFMT -using Random: Sampler, SamplerRangeFast, SamplerRangeInt, MT_CACHE_F, MT_CACHE_I +using Random: Sampler, SamplerRangeFast, SamplerRangeInt, SamplerRangeNDL, MT_CACHE_F, MT_CACHE_I import Future # randjump @@ -222,25 +222,11 @@ randmtzig_fill_ziggurat_tables() @test all(we == Random.we) @test all(fe == Random.fe) -#same random numbers on for small ranges on all systems -guardseed() do - seed = rand(UInt) - Random.seed!(seed) - r = map(Int64, rand(map(Int32, 97:122))) - Random.seed!(seed) - @test r == rand(map(Int64, 97:122)) - Random.seed!(seed) - r = map(UInt64, rand(map(UInt32, 97:122))) - Random.seed!(seed) - @test r == rand(map(UInt64, 97:122)) -end - for U in (Int64, UInt64) @test all(div(one(UInt64) << 52, k)*k - 1 == SamplerRangeInt(map(U, 1:k)).u for k in 13 .+ Int64(2).^(1:30)) end - #issue 8257 let i8257 = 1:1/3:100 for i = 1:100 @@ -707,11 +693,16 @@ struct RandomStruct23964 end @test_throws ArgumentError rand(RandomStruct23964()) end -@testset "rand(::$(typeof(RNG)), ::UnitRange{$T}" for RNG ∈ (MersenneTwister(), RandomDevice()), - T ∈ (Int32, UInt32, Int64, Int128, UInt128) - RNG isa MersenneTwister && Random.seed!(RNG, rand(UInt128)) # for reproducibility - r = T(1):T(108) - @test rand(RNG, SamplerRangeFast(r)) ∈ r +@testset "rand(::$(typeof(RNG)), ::UnitRange{$T}" for RNG ∈ (MersenneTwister(rand(UInt128)), RandomDevice()), + T ∈ (Int8, Int16, Int32, UInt32, Int64, Int128, UInt128) + for S in (SamplerRangeInt, SamplerRangeFast, SamplerRangeNDL) + S == SamplerRangeNDL && sizeof(T) > 8 && continue + r = T(1):T(108) + @test rand(RNG, S(r)) ∈ r + @test rand(RNG, S(typemin(T):typemax(T))) isa T + a, b = sort!(rand(-1000:1000, 2) .% T) + @test rand(RNG, S(a:b)) ∈ a:b + end end @testset "rand! is allocation-free" begin @@ -803,8 +794,8 @@ end @test rand!(GLOBAL_RNG, A, x) === A == rand!(mt, B, x) === B end # issue #33170 - @test Sampler(GLOBAL_RNG, 2:4, Val(1)) isa SamplerRangeFast - @test Sampler(GLOBAL_RNG, 2:4, Val(Inf)) isa SamplerRangeFast + @test Sampler(GLOBAL_RNG, 2:4, Val(1)) isa SamplerRangeNDL + @test Sampler(GLOBAL_RNG, 2:4, Val(Inf)) isa SamplerRangeNDL end @testset "RNGs broadcast as scalars: T" for T in (MersenneTwister, RandomDevice)