implement "nearly division less" algorithm for rand(a: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
 
 

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

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
 ```

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.

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
 

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

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)