From 5facdd1feb6c0122af6cd26b17fbecc45d835fed Mon Sep 17 00:00:00 2001
From: Rafael Fourquet <fourquet.rafael@gmail.com>
Date: Tue, 18 Sep 2018 11:52:05 +0200
Subject: [PATCH 1/7] implement "nearly division less" algorithm for rand(a:b)

Cf. https://arxiv.org/abs/1805.10941.
Closes #20582, #29004.
---
 stdlib/Random/src/generation.jl | 51 +++++++++++++++++++++++++++++----
 stdlib/Random/test/runtests.jl  |  9 ++++++
 2 files changed, 55 insertions(+), 5 deletions(-)

diff --git a/stdlib/Random/src/generation.jl b/stdlib/Random/src/generation.jl
index caa5e3cde3f45..d1e03266fc162 100644
--- a/stdlib/Random/src/generation.jl
+++ b/stdlib/Random/src/generation.jl
@@ -166,13 +166,15 @@ 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
+# 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", 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) "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
+# 3) "Nearly Division Less", which is the fastest algorithm for types of size
+#    up to 64 bits. This will be the default in a future realease.
 
 #### helper functions
 
@@ -295,6 +297,45 @@ 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
+
+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)
+        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
+
+# API until this is made the default
+fast(r::AbstractUnitRange{<:Base.BitInteger64}) = SamplerRangeNDL(r)
+fast(r::AbstractArray) = Random.SamplerSimple(r, fast(firstindex(r):lastindex(r)))
+
 
 ### BigInt
 
diff --git a/stdlib/Random/test/runtests.jl b/stdlib/Random/test/runtests.jl
index 20f81f73e344a..1db2350702b3b 100644
--- a/stdlib/Random/test/runtests.jl
+++ b/stdlib/Random/test/runtests.jl
@@ -810,3 +810,12 @@ end
 @testset "RNGs broadcast as scalars: T" for T in (MersenneTwister, RandomDevice)
     @test length.(rand.(T(), 1:3)) == 1:3
 end
+
+@testset "fast(a:b)" begin
+    for bounds = (rand(Int, 2), rand(-1000:1000, 2))
+        a, b = minmax(bounds...)
+        a == b && continue
+        @test rand(Random.fast(a:b)) ∈ a:b
+        @test rand(Random.fast(a:1:b)) ∈ a:b
+    end
+end

From 60c3edd5f950898357e1141c8b98e676348936af Mon Sep 17 00:00:00 2001
From: Rafael Fourquet <fourquet.rafael@gmail.com>
Date: Tue, 18 Sep 2018 15:53:45 +0200
Subject: [PATCH 2/7] fix overflow error in tests

---
 stdlib/Random/test/runtests.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/stdlib/Random/test/runtests.jl b/stdlib/Random/test/runtests.jl
index 1db2350702b3b..7aaeed660bb74 100644
--- a/stdlib/Random/test/runtests.jl
+++ b/stdlib/Random/test/runtests.jl
@@ -813,7 +813,7 @@ end
 
 @testset "fast(a:b)" begin
     for bounds = (rand(Int, 2), rand(-1000:1000, 2))
-        a, b = minmax(bounds...)
+        a, b = minmax((bounds .>> 2)...) # right-shift to avoid overflow in length(a:1:b)
         a == b && continue
         @test rand(Random.fast(a:b)) ∈ a:b
         @test rand(Random.fast(a:1:b)) ∈ a:b

From e4ea174d51c85626cbdc4dc1f7f055a0dd8334e3 Mon Sep 17 00:00:00 2001
From: Rafael Fourquet <fourquet.rafael@gmail.com>
Date: Fri, 31 May 2019 11:44:16 +0200
Subject: [PATCH 3/7] make NDL the default algo

---
 stdlib/Random/src/RNGs.jl       |  8 ------
 stdlib/Random/src/generation.jl | 20 +++++++++------
 stdlib/Random/test/runtests.jl  | 44 ++++++++++-----------------------
 3 files changed, 25 insertions(+), 47 deletions(-)

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 d1e03266fc162..717b24ac9921b 100644
--- a/stdlib/Random/src/generation.jl
+++ b/stdlib/Random/src/generation.jl
@@ -171,10 +171,18 @@ end
 # 1) "Fast", 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) "Default" which tries to use as few entropy bits as possible, at the cost of a
+# 2) "Original" 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
-# 3) "Nearly Division Less", which is the fastest algorithm for types of size
-#    up to 64 bits. This will be the default in a future realease.
+# 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
 
@@ -226,7 +234,7 @@ function rand(rng::AbstractRNG, sp::SamplerRangeFast{UInt128,T}) where T
     x % T + a
 end
 
-#### Default
+#### Original
 
 # remainder function according to Knuth, where rem_knuth(a, 0) = a
 rem_knuth(a::UInt, b::UInt) = a % (b + (b == 0)) + a * (b == 0)
@@ -276,10 +284,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
 
diff --git a/stdlib/Random/test/runtests.jl b/stdlib/Random/test/runtests.jl
index 7aaeed660bb74..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,19 +794,10 @@ 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)
     @test length.(rand.(T(), 1:3)) == 1:3
 end
-
-@testset "fast(a:b)" begin
-    for bounds = (rand(Int, 2), rand(-1000:1000, 2))
-        a, b = minmax((bounds .>> 2)...) # right-shift to avoid overflow in length(a:1:b)
-        a == b && continue
-        @test rand(Random.fast(a:b)) ∈ a:b
-        @test rand(Random.fast(a:1:b)) ∈ a:b
-    end
-end

From bb402fe4972e6f5870c5ce312d24277f8db0d3ca Mon Sep 17 00:00:00 2001
From: Rafael Fourquet <fourquet.rafael@gmail.com>
Date: Sat, 11 Apr 2020 17:31:16 +0200
Subject: [PATCH 4/7] update NEWS.md

---
 NEWS.md | 7 +++++++
 1 file changed, 7 insertions(+)

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
 
 

From f291355b5c580d5f3edad13adef3800119767c7a Mon Sep 17 00:00:00 2001
From: Rafael Fourquet <fourquet.rafael@gmail.com>
Date: Mon, 13 Apr 2020 09:29:15 +0200
Subject: [PATCH 5/7] try make tests pass on 32-bits machines

---
 stdlib/LinearAlgebra/test/bunchkaufman.jl |  2 +-
 stdlib/Random/docs/src/index.md           | 18 +++++++++---------
 stdlib/Random/src/misc.jl                 | 10 +++++-----
 3 files changed, 15 insertions(+), 15 deletions(-)

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

From fab8970af0d0dad1bba3829d1780cff20b5bc074 Mon Sep 17 00:00:00 2001
From: Rafael Fourquet <fourquet.rafael@gmail.com>
Date: Mon, 27 Apr 2020 15:34:44 +0200
Subject: [PATCH 6/7] add a comment for mod(-s, s)

---
 stdlib/Random/src/generation.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/stdlib/Random/src/generation.jl b/stdlib/Random/src/generation.jl
index 717b24ac9921b..ee1d8db5c62c9 100644
--- a/stdlib/Random/src/generation.jl
+++ b/stdlib/Random/src/generation.jl
@@ -303,7 +303,7 @@ end
 
 #### Nearly Division Less
 
-# cf. https://arxiv.org/abs/1805.10941
+# cf. https://arxiv.org/abs/1805.10941 (algorithm 5)
 
 struct SamplerRangeNDL{U<:Unsigned,T} <: Sampler{T}
     a::T  # first element of the range
@@ -326,7 +326,7 @@ function rand(rng::AbstractRNG, sp::SamplerRangeNDL{U,T}) where {U,T}
     m = x * s
     l = m % U
     if l < s
-        t = mod(-s, 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

From b7348d427e1de6b7ea7d33bbd3e4599728ca4e3b Mon Sep 17 00:00:00 2001
From: Rafael Fourquet <fourquet.rafael@gmail.com>
Date: Fri, 1 May 2020 12:17:58 +0200
Subject: [PATCH 7/7] remove vestigial transient `fast` function, and update
 comments

---
 stdlib/Random/src/generation.jl | 17 +++++++----------
 1 file changed, 7 insertions(+), 10 deletions(-)

diff --git a/stdlib/Random/src/generation.jl b/stdlib/Random/src/generation.jl
index ee1d8db5c62c9..329620c3ab40b 100644
--- a/stdlib/Random/src/generation.jl
+++ b/stdlib/Random/src/generation.jl
@@ -168,11 +168,12 @@ end
 
 # 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", 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) "Original" 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
+# 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
@@ -234,7 +235,7 @@ function rand(rng::AbstractRNG, sp::SamplerRangeFast{UInt128,T}) where T
     x % T + a
 end
 
-#### Original
+#### "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)
@@ -336,10 +337,6 @@ function rand(rng::AbstractRNG, sp::SamplerRangeNDL{U,T}) where {U,T}
     (s == 0 ? x : m >> (8*sizeof(U))) % T + sp.a
 end
 
-# API until this is made the default
-fast(r::AbstractUnitRange{<:Base.BitInteger64}) = SamplerRangeNDL(r)
-fast(r::AbstractArray) = Random.SamplerSimple(r, fast(firstindex(r):lastindex(r)))
-
 
 ### BigInt