make TaskLocal the default RNG

JeffBezanson · JeffBezanson · commit 0e3ea6e4352a · 2021-06-01T15:16:06.000-04:00
diff --git a/NEWS.md b/NEWS.md
@@ -26,6 +26,8 @@ Language changes
 * Destructuring will no longer mutate values on the left hand side while iterating through values on the right hand side. In the example
   of an array `x`, `x[2], x[1] = x` will now swap the first and second entry of `x`, whereas it used to fill both entries with `x[1]`
   because `x[2]` was mutated during the iteration of `x`. ([#40737])
+* The default random number generator has changed, so all random numbers will be different (even with the
+  same seed) unless an explicit RNG object is used. See the section on the `Random` standard library below.
 
 Compiler/Runtime improvements
 -----------------------------
@@ -39,7 +41,12 @@ Command-line option changes
 
 Multi-threading changes
 -----------------------
+
 * If the `JULIA_NUM_THREADS` environment variable is set to `auto`, then the number of threads will be set to the number of CPU threads ([#38952])
+* Every `Task` object has a local random number generator state, providing reproducible (schedule-independent) execution
+  of parallel simulation code by default. The default generator is also significantly faster in parallel than in
+  previous versions.
+
 
 Build system changes
 --------------------
@@ -130,6 +137,9 @@ Standard library changes
 
 #### Random
 
+* The default random number generator has been changed from Mersenne Twister to [Xoshiro256++](https://prng.di.unimi.it/).
+  The new generator has smaller state, better performance, and superior statistical properties.
+  This generator is the one used for reproducible Task-local randomness.
 
 #### REPL
 
diff --git a/doc/src/devdocs/subarrays.md b/doc/src/devdocs/subarrays.md
@@ -19,14 +19,14 @@ julia> A = rand(2,3,4);
 
 julia> S1 = view(A, :, 1, 2:3)
 2×2 view(::Array{Float64, 3}, :, 1, 2:3) with eltype Float64:
- 0.200586  0.066423
- 0.298614  0.956753
+ 0.98627   0.991588
+ 0.523355  0.105192
 
 julia> S2 = view(A, 1, :, 2:3)
 3×2 view(::Array{Float64, 3}, 1, :, 2:3) with eltype Float64:
- 0.200586  0.066423
- 0.246837  0.646691
- 0.648882  0.276021
+ 0.98627    0.991588
+ 0.0890391  0.904316
+ 0.134982   0.299203
 ```
 ```@meta
 DocTestSetup = nothing
diff --git a/doc/src/manual/performance-tips.md b/doc/src/manual/performance-tips.md
@@ -77,12 +77,12 @@ julia> function sum_global()
        end;
 
 julia> @time sum_global()
-  0.009639 seconds (7.36 k allocations: 300.310 KiB, 98.32% compilation time)
-496.84883432553846
+  0.026328 seconds (9.30 k allocations: 416.747 KiB, 36.50% gc time, 99.48% compilation time)
+508.39048990953665
 
 julia> @time sum_global()
-  0.000140 seconds (3.49 k allocations: 70.313 KiB)
-496.84883432553846
+  0.000075 seconds (3.49 k allocations: 70.156 KiB)
+508.39048990953665
 ```
 
 On the first call (`@time sum_global()`) the function gets compiled. (If you've not yet used [`@time`](@ref)
@@ -113,12 +113,12 @@ julia> function sum_arg(x)
        end;
 
 julia> @time sum_arg(x)
-  0.006202 seconds (4.18 k allocations: 217.860 KiB, 99.72% compilation time)
-496.84883432553846
+  0.010298 seconds (4.23 k allocations: 226.021 KiB, 99.81% compilation time)
+508.39048990953665
 
 julia> @time sum_arg(x)
   0.000005 seconds (1 allocation: 16 bytes)
-496.84883432553846
+508.39048990953665
 ```
 
 The 1 allocation seen is from running the `@time` macro itself in global scope. If we instead run
@@ -129,7 +129,7 @@ julia> time_sum(x) = @time sum_arg(x);
 
 julia> time_sum(x)
   0.000001 seconds
-496.84883432553846
+508.39048990953665
 ```
 
 In some situations, your function may need to allocate memory as part of its operation, and this
diff --git 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!(12343210)
+Random.seed!(12343212)
 
 areal = randn(n,n)/2
 aimg  = randn(n,n)/2
diff --git a/stdlib/LinearAlgebra/test/cholesky.jl b/stdlib/LinearAlgebra/test/cholesky.jl
@@ -39,7 +39,7 @@ end
     n1 = div(n, 2)
     n2 = 2*n1
 
-    Random.seed!(12343)
+    Random.seed!(12344)
 
     areal = randn(n,n)/2
     aimg  = randn(n,n)/2
diff --git a/stdlib/LinearAlgebra/test/dense.jl b/stdlib/LinearAlgebra/test/dense.jl
@@ -15,17 +15,17 @@ n = 10
 n1 = div(n, 2)
 n2 = 2*n1
 
-Random.seed!(1234321)
+Random.seed!(1234323)
 
 @testset "Matrix condition number" begin
     ainit = rand(n,n)
     @testset "for $elty" for elty in (Float32, Float64, ComplexF32, ComplexF64)
         ainit = convert(Matrix{elty}, ainit)
         for a in (copy(ainit), view(ainit, 1:n, 1:n))
-            @test cond(a,1) ≈ 4.837320054554436e+02 atol=0.01
-            @test cond(a,2) ≈ 1.960057871514615e+02 atol=0.01
-            @test cond(a,Inf) ≈ 3.757017682707787e+02 atol=0.01
-            @test cond(a[:,1:5]) ≈ 10.233059337453463 atol=0.01
+            @test cond(a,1) ≈ 162.93076496229205 atol=0.5
+            @test cond(a,2) ≈ 97.7362909028399 atol=0.5
+            @test cond(a,Inf) ≈ 196.00888526411242 atol=0.4
+            @test cond(a[:,1:5]) ≈ 10.466183798451308 atol=0.01
             @test_throws ArgumentError cond(a,3)
         end
     end
diff --git a/stdlib/LinearAlgebra/test/eigen.jl b/stdlib/LinearAlgebra/test/eigen.jl
@@ -11,7 +11,7 @@ n = 10
 n1 = div(n, 2)
 n2 = 2*n1
 
-Random.seed!(1234321)
+Random.seed!(12343215)
 
 areal = randn(n,n)/2
 aimg  = randn(n,n)/2
diff --git a/stdlib/LinearAlgebra/test/lu.jl b/stdlib/LinearAlgebra/test/lu.jl
@@ -11,7 +11,7 @@ n = 10
 n1 = div(n, 2)
 n2 = 2*n1
 
-Random.seed!(1234321)
+Random.seed!(1234323)
 
 areal = randn(n,n)/2
 aimg  = randn(n,n)/2
diff --git a/stdlib/LinearAlgebra/test/qr.jl b/stdlib/LinearAlgebra/test/qr.jl
@@ -11,7 +11,7 @@ n = 10
 n1 = div(n, 2)
 n2 = 2*n1
 
-Random.seed!(1234321)
+Random.seed!(1234325)
 
 areal = randn(n,n)/2
 aimg  = randn(n,n)/2
diff --git a/stdlib/LinearAlgebra/test/uniformscaling.jl b/stdlib/LinearAlgebra/test/uniformscaling.jl
@@ -10,7 +10,7 @@ using .Main.Quaternions
 isdefined(Main, :OffsetArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "OffsetArrays.jl"))
 using .Main.OffsetArrays
 
-Random.seed!(123)
+Random.seed!(1234543)
 
 @testset "basic functions" begin
     @test I === I' # transpose
diff --git a/stdlib/Random/docs/src/index.md b/stdlib/Random/docs/src/index.md
@@ -147,20 +147,20 @@ Scalar and array methods for `Die` now work as expected:
 
 ```jldoctest Die; setup = :(Random.seed!(1))
 julia> rand(Die)
-Die(15)
+Die(6)
 
 julia> rand(MersenneTwister(0), Die)
 Die(11)
 
 julia> rand(Die, 3)
 3-element Vector{Die}:
- Die(18)
- Die(5)
+ Die(15)
+ Die(19)
  Die(4)
 
 julia> a = Vector{Die}(undef, 3); rand!(a)
 3-element Vector{Die}:
- Die(5)
+ Die(17)
  Die(20)
  Die(15)
 ```
@@ -175,13 +175,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))
-3
+1
 
 julia> rand(Die(4), 3)
 3-element Vector{Any}:
+ 3
  4
  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
@@ -355,48 +355,37 @@ function seed!(r::MersenneTwister, seed::Vector{UInt32})
     return r
 end
 
-seed!(r::MersenneTwister=default_rng()) = seed!(r, make_seed())
+seed!(r::MersenneTwister) = seed!(r, make_seed())
 seed!(r::MersenneTwister, n::Integer) = seed!(r, make_seed(n))
-seed!(seed::Union{Integer,Vector{UInt32}}) = seed!(default_rng(), seed)
 
 
 ### Global RNG
 
-const THREAD_RNGs = MersenneTwister[]
-@inline default_rng() = default_rng(Threads.threadid())
-@noinline function default_rng(tid::Int)
-    0 < tid <= length(THREAD_RNGs) || _rng_length_assert()
-    if @inbounds isassigned(THREAD_RNGs, tid)
-        @inbounds MT = THREAD_RNGs[tid]
-    else
-        MT = MersenneTwister()
-        @inbounds THREAD_RNGs[tid] = MT
-    end
-    return MT
-end
-@noinline _rng_length_assert() =  @assert false "0 < tid <= length(THREAD_RNGs)"
-
-function __init__()
-    resize!(empty!(THREAD_RNGs), Threads.nthreads()) # ensures that we didn't save a bad object
-    seed!(TaskLocalRNG())
-end
-
-
 struct _GLOBAL_RNG <: AbstractRNG
     global const GLOBAL_RNG = _GLOBAL_RNG.instance
 end
 
-# GLOBAL_RNG currently represents a MersenneTwister
-typeof_rng(::_GLOBAL_RNG) = MersenneTwister
+# GLOBAL_RNG currently uses TaskLocalRNG
+typeof_rng(::_GLOBAL_RNG) = TaskLocalRNG
 
-copy!(dst::MersenneTwister, ::_GLOBAL_RNG) = copy!(dst, default_rng())
-copy!(::_GLOBAL_RNG, src::MersenneTwister) = copy!(default_rng(), src)
+@inline default_rng() = TaskLocalRNG()
+@inline default_rng(tid::Int) = TaskLocalRNG()
+
+copy!(dst::Xoshiro, ::_GLOBAL_RNG) = copy!(dst, default_rng())
+copy!(::_GLOBAL_RNG, src::Xoshiro) = copy!(default_rng(), src)
 copy(::_GLOBAL_RNG) = copy(default_rng())
 
-seed!(::_GLOBAL_RNG, seed::Vector{UInt32}) = seed!(default_rng(), seed)
-seed!(::_GLOBAL_RNG, n::Integer) = seed!(default_rng(), n)
-seed!(::_GLOBAL_RNG, ::Nothing) = seed!(default_rng(), nothing)
-seed!(::_GLOBAL_RNG) = seed!(default_rng(), nothing)
+GLOBAL_SEED = 0
+
+seed!(::_GLOBAL_RNG, seed) = (global GLOBAL_SEED = seed; seed!(default_rng(), seed))
+
+function seed!(rng::_GLOBAL_RNG)
+    seed!(rng, (rand(RandomDevice(), UInt64), rand(RandomDevice(), UInt64),
+                rand(RandomDevice(), UInt64), rand(RandomDevice(), UInt64)))
+end
+seed!() = seed!(GLOBAL_RNG)
+seed!(rng::_GLOBAL_RNG, ::Nothing) = seed!(rng)  # to resolve ambiguity
+seed!(seed::Union{Integer,Vector{UInt32},Vector{UInt64},NTuple{4,UInt64}}) = seed!(GLOBAL_RNG, seed)
 
 rng_native_52(::_GLOBAL_RNG) = rng_native_52(default_rng())
 rand(::_GLOBAL_RNG, sp::SamplerBoolBitInteger) = rand(default_rng(), sp)
@@ -420,6 +409,10 @@ for T in BitInteger_types
     @eval rand!(::_GLOBAL_RNG, A::Array{$T}, I::SamplerType{$T}) = rand!(default_rng(), A, I)
 end
 
+function __init__()
+    seed!(GLOBAL_RNG)
+end
+
 
 ### generation
 
diff --git a/stdlib/Random/src/misc.jl b/stdlib/Random/src/misc.jl
@@ -53,13 +53,13 @@ number generator, see [Random Numbers](@ref).
 # Examples
 ```jldoctest
 julia> Random.seed!(3); randstring()
-"Y7m62wOj"
+"vZmAMp3z"
 
 julia> randstring(MersenneTwister(3), 'a':'z', 6)
 "ocucay"
 
 julia> randstring("ACGT")
-"ATTTGCGT"
+"CAAACACC"
 ```
 
 !!! note
diff --git a/stdlib/Random/test/runtests.jl b/stdlib/Random/test/runtests.jl
@@ -628,7 +628,7 @@ guardseed() do
     m = MersenneTwister(0)
     @test Random.seed!() === g
     @test Random.seed!(rand(UInt)) === g
-    @test Random.seed!(rand(UInt32, rand(1:10))) === g
+    @test Random.seed!(rand(UInt32, rand(1:8))) === g
     @test Random.seed!(m) === m
     @test Random.seed!(m, rand(UInt)) === m
     @test Random.seed!(m, rand(UInt32, rand(1:10))) === m
@@ -752,28 +752,26 @@ end
     @test Random.seed!(GLOBAL_RNG, 0) === LOCAL_RNG
     @test Random.seed!(GLOBAL_RNG) === LOCAL_RNG
 
-    mt = MersenneTwister(1)
-    @test copy!(mt, GLOBAL_RNG) === mt
-    @test mt == LOCAL_RNG
-    Random.seed!(mt, 2)
-    @test mt != LOCAL_RNG
-    @test copy!(GLOBAL_RNG, mt) === LOCAL_RNG
-    @test mt == LOCAL_RNG
-    mt2 = copy(GLOBAL_RNG)
-    @test mt2 isa typeof(LOCAL_RNG)
-    @test mt2 !== LOCAL_RNG
-    @test mt2 == LOCAL_RNG
+    xo = Xoshiro()
+    @test copy!(xo, GLOBAL_RNG) === xo
+    @test xo == LOCAL_RNG
+    Random.seed!(xo, 2)
+    @test xo != LOCAL_RNG
+    @test copy!(GLOBAL_RNG, xo) === LOCAL_RNG
+    @test xo == LOCAL_RNG
+    xo2 = copy(GLOBAL_RNG)
+    @test xo2 !== LOCAL_RNG
+    @test xo2 == LOCAL_RNG
 
     for T in (Random.UInt52Raw{UInt64},
-              Random.UInt2x52Raw{UInt128},
               Random.UInt104Raw{UInt128},
               Random.CloseOpen12_64)
         x = Random.SamplerTrivial(T())
-        @test rand(GLOBAL_RNG, x) === rand(mt, x)
+        @test rand(GLOBAL_RNG, x) === rand(xo, x)
     end
     for T in (Int64, UInt64, Int128, UInt128, Bool, Int8, UInt8, Int16, UInt16, Int32, UInt32)
         x = Random.SamplerType{T}()
-        @test rand(GLOBAL_RNG, x) === rand(mt, x)
+        @test rand(GLOBAL_RNG, x) === rand(xo, x)
     end
 
     A = fill(0.0, 100, 100)
@@ -782,21 +780,21 @@ end
     vB = view(B, :, :)
     I1 = Random.SamplerTrivial(Random.CloseOpen01{Float64}())
     I2 = Random.SamplerTrivial(Random.CloseOpen12{Float64}())
-    @test rand!(GLOBAL_RNG, A, I1) === A == rand!(mt, B, I1) === B
+    @test rand!(GLOBAL_RNG, A, I1) === A == rand!(xo, B, I1) === B
     B = fill!(B, 1.0)
     @test rand!(GLOBAL_RNG, vA, I1) === vA
-    rand!(mt, vB, I1)
+    rand!(xo, vB, I1)
     @test A == B
     for T in (Float16, Float32)
         B = fill!(B, 1.0)
-        @test rand!(GLOBAL_RNG, A, I2) === A == rand!(mt, B, I2) === B
+        @test rand!(GLOBAL_RNG, A, I2) === A == rand!(xo, B, I2) === B
         B = fill!(B, 1.0)
-        @test rand!(GLOBAL_RNG, A, I1) === A == rand!(mt, B, I1) === B
+        @test rand!(GLOBAL_RNG, A, I1) === A == rand!(xo, B, I1) === B
     end
     for T in Base.BitInteger_types
         x = Random.SamplerType{T}()
         B = fill!(B, 1.0)
-        @test rand!(GLOBAL_RNG, A, x) === A == rand!(mt, B, x) === B
+        @test rand!(GLOBAL_RNG, A, x) === A == rand!(xo, B, x) === B
     end
     # issue #33170
     @test Sampler(GLOBAL_RNG, 2:4, Val(1)) isa SamplerRangeNDL
diff --git a/stdlib/SparseArrays/src/sparsematrix.jl b/stdlib/SparseArrays/src/sparsematrix.jl
diff --git a/stdlib/SparseArrays/test/sparse.jl b/stdlib/SparseArrays/test/sparse.jl
diff --git a/stdlib/Test/src/Test.jl b/stdlib/Test/src/Test.jl
diff --git a/stdlib/Test/test/runtests.jl b/stdlib/Test/test/runtests.jl
diff --git a/test/error.jl b/test/error.jl
