ditch fork!, add array version, do not support TaskLocalRNG

Now that seed! accepts an rng as 2nd arg, fork! is replaced by seed!.
Array version is useful for `@threads` (cf. docstring).
`fork(TaskLocalRNG())` is ambiguous: return a Xoshiro or TaskLocalRNG() ?
Leave that question for later.
Also, returning Xoshiro was seriously
complicating the docstring.

Author: rfourquet

NEWS.md

@@ -49,7 +49,7 @@ Standard library changes
 
 #### Random
 
-* It's now possible to efficiently create a new `Xoshiro` RNG from an existing one via `Random.fork`, which
+* It's now possible to efficiently create a new `Xoshiro` instance from an existing one via `Random.fork`, which
   can be useful in parallel computations ([#58193]).
 
 #### REPL

stdlib/Random/docs/src/index.md

@@ -83,7 +83,6 @@ Random.TaskLocalRNG
 Random.Xoshiro
 Random.MersenneTwister
 Random.RandomDevice
-Random.fork!
 Random.fork
 ```
 

stdlib/Random/src/RNGs.jl

@@ -306,3 +306,61 @@ function hash_seed(str::AbstractString, ctx::SHA_CTX)
     end
     SHA.update!(ctx, (0x05,))
 end
+
+
+## forking
+
+"""
+    Random.fork(rng::AbstractRNG, [dims...])::typeof(rng)
+
+Create deterministically a new independent RNG object from an existing one, of the same type.
+When `dims` is specified (as integers or a tuple of integers),
+an array of such independent RNGs is created.
+
+This is the recommended way to initialize fresh RNG instances when reproducibility might be
+required, and can be useful in particular in multi-threaded contexts, where race conditions
+must be avoided. For example:
+```julia
+function dotask(rng::Xoshiro)
+    num_subtasks = 10
+    rngs = Random.fork(rng, num_subtasks)
+    result = Vector(undef, num_subtasks)
+    Threads.@threads for i=1:num_subtasks
+        result[i] = dosubtask(i, rngs[i])
+    end
+    result
+end
+```
+Note that this functions generally mutates its `rng` argument, so calling it on the same `rng`
+must not be done concurrently from multiple threads.
+
+Currently, a method is only implemented for `Xoshiro`.
+
+!!! note
+    `Random.fork(::Xoshiro)` uses the same algorithm as when the task-local RNG
+    of a new task is created from the task-local RNG of the parent task.
+
+!!! compat "Julia 1.13"
+    This function was introduced in Julia 1.13.
+
+# Examples
+```jldoctest
+julia> x = Xoshiro(0)
+Xoshiro(0xdb2fa90498613fdf, 0x48d73dc42d195740, 0x8c49bc52dc8a77ea, 0x1911b814c02405e8, 0x22a21880af5dc689)
+
+julia> [x, fork(x)] # x is mutated
+2-element Vector{Xoshiro}:
+ Xoshiro(0xdb2fa90498613fdf, 0x48d73dc42d195740, 0x8c49bc52dc8a77ea, 0x1911b814c02405e8, 0x26b589433d8074be)
+ Xoshiro(0x545f53b997598dfb, 0xac80f92d91bb35a5, 0xb6eb3382e8c409ca, 0xa9aa6968fbdd5e83, 0x26b589433d8074be)
+
+julia> fork(x, 3)
+3-element Vector{Xoshiro}:
+ Xoshiro(0xfc86733bafa6df6d, 0x5ff051ecb5937fcf, 0x935c4e55a82ca686, 0xa57b44768cdb84e9, 0x845f4ebfc53d5497)
+ Xoshiro(0x9c49327a59542654, 0x7c2f821b7716e6b7, 0x586a3fe58fed92f7, 0x28bbf526c1aca281, 0x425e2dc6f55934e4)
+ Xoshiro(0xd5cbe598083243e0, 0xfba09a94aaf998af, 0xa674c207f3796c54, 0x084d4986ad49c4eb, 0x52a722b1a914a4b5)
+```
+"""
+fork
+
+fork(rng::AbstractRNG, dims::Dims) = typeof(rng)[fork(rng) for _=LinearIndices(dims)]
+fork(rng::AbstractRNG, d1::Integer, dims::Integer...) = fork(rng, Dims((d1, dims...)))

stdlib/Random/src/Random.jl

@@ -29,7 +29,7 @@ export rand!, randn!,
        randcycle, randcycle!,
        AbstractRNG, MersenneTwister, RandomDevice, TaskLocalRNG, Xoshiro
 
-public fork, fork!, seed!, default_rng, Sampler, SamplerType, SamplerTrivial, SamplerSimple
+public fork, seed!, default_rng, Sampler, SamplerType, SamplerTrivial, SamplerSimple
 
 ## general definitions
 

stdlib/Random/src/Xoshiro.jl

@@ -246,6 +246,9 @@ hash(x::Union{TaskLocalRNG, Xoshiro}, h::UInt) = hash(getstate(x), h + 0x49a62c2
 seed!(rng::Union{TaskLocalRNG, Xoshiro}, seeder::AbstractRNG) =
     initstate!(rng, rand(seeder, NTuple{4, UInt64}))
 
+# when seeder is a Xoshiro, use forking instead of regular seeding, to avoid correlations
+seed!(rng::Union{TaskLocalRNG, Xoshiro}, seeder::Union{TaskLocalRNG, Xoshiro}) =
+    setstate!(rng, _fork(seeder))
 
 @inline function rand(x::Union{TaskLocalRNG, Xoshiro}, ::SamplerType{UInt64})
     s0, s1, s2, s3 = getstate(x)
@@ -332,33 +335,4 @@ function _fork(src::Union{Xoshiro, TaskLocalRNG})
     (state..., s4)
 end
 
-"""
-    Random.fork!(dst::Union{Xoshiro, TaskLocalRNG}, src::Union{Xoshiro, TaskLocalRNG} = TaskLocalRNG()) -> dst
-
-Equivalent to `copy!(dst, fork(src))`.
-See also [`fork`](@ref).
-
-!!! compat "Julia 1.13"
-    This function was introduced in Julia 1.13.
-"""
-fork!(dst::Union{Xoshiro, TaskLocalRNG}, src::Union{Xoshiro, TaskLocalRNG}=TaskLocalRNG()) =
-    setstate!(dst, _fork(src))
-
-"""
-    Random.fork(src::Union{Xoshiro, TaskLocalRNG} = TaskLocalRNG())::Xoshiro
-
-Create a new `Xoshiro` object from `src`, in the same way that the task local RNG of a new
-task is created from the task local RNG of the parent task.
-This is the recommended way to initialize a fresh RNG from an existing one.
-
-!!! note
-    When `src` is of type `TaskLocalRNG`, this function is guaranteed to return an RNG of
-    type `Xoshiro` only as long as `Xoshiro` and `TaskLocalRNG` are implementing the same
-    underlying generator. It may change in the future.
-
-!!! compat "Julia 1.13"
-    This function was introduced in Julia 1.13.
-
-See also [`Random.fork!`](@ref).
-"""
-fork(src::Union{Xoshiro, TaskLocalRNG}=TaskLocalRNG()) = Xoshiro(_fork(src)...)
+fork(src::Xoshiro) = Xoshiro(_fork(src)...)

stdlib/Random/test/runtests.jl

@@ -1302,40 +1302,57 @@ end
 
 @testset "fork" begin
     xx = copy(TaskLocalRNG())
+    x0 = copy(xx)
     x1 = Random.fork(xx)
     x2 = fetch(@async copy(TaskLocalRNG()))
     @test x1 isa Xoshiro && x2 isa Xoshiro
     @test x1 == x2 # currently, equality involves all 5 UInt64 words of the state
     @test xx == TaskLocalRNG()
 
-    x3 = Random.fork(TaskLocalRNG())
-    @test x3 isa Xoshiro
-    copy!(TaskLocalRNG(), xx) # reset its state
-    x4 = Random.fork(xx)
-    @test x4 isa Xoshiro
-    @test x3 == x4
-    copy!(xx, TaskLocalRNG())
-
-    @test xx == TaskLocalRNG() # check assumptions
-    x5 = Random.fork()
-    @test xx != TaskLocalRNG() # TaskLocalRNG() was forked off
-    copy!(TaskLocalRNG(), xx)
-    @test x5 == x4
-
-    x6 = Xoshiro(0, 0, 0, 0, 0)
-    @test x6 === Random.fork!(x6, xx)
-    copy!(xx, TaskLocalRNG())
-    @test x6 == x5
-    @test x6 === Random.fork!(x6, TaskLocalRNG())
-    copy!(TaskLocalRNG(), xx)
-    @test x6 == x5
-    @test xx == TaskLocalRNG() # check assumptions
-    @test x6 === Random.fork!(x6)
-    @test xx != TaskLocalRNG()
-    copy!(TaskLocalRNG(), xx)
-    @test x6 == x5
-
-    @test TaskLocalRNG() === Random.fork!(TaskLocalRNG(), copy(xx))
-    @test x6 === Random.fork!(x6, copy(xx))
-    @test TaskLocalRNG() == x6
+    # fork is deterministic
+    copy!(xx, x0)
+    x3 = Random.fork(xx)
+    @test x3 == x2
+    @test rand(x3, UInt128) == rand(x2, UInt128)
+
+    # seed! uses the same mechanism
+    copy!(xx, x0)
+    x4 = Xoshiro(0, 0, 0, 0, 0)
+    @test x4 === Random.seed!(x4, xx)
+    @test x4 == x1
+    copy!(TaskLocalRNG(), x0)
+    x5 = Xoshiro(0, 0, 0, 0, 0)
+    @test x5 === Random.seed!(x5, TaskLocalRNG())
+    @test x5 == x1
+    @test xx == TaskLocalRNG() # both are in the same state after being forked off
+
+    @test TaskLocalRNG() == Random.seed!(TaskLocalRNG(), copy!(xx, x0))
+    @test TaskLocalRNG() == x4
+    copy!(TaskLocalRNG(), x0)
+    @test TaskLocalRNG() == Random.seed!(TaskLocalRNG(), TaskLocalRNG())
+
+    # self-seeding
+    copy!(xx, x0)
+    copy!(TaskLocalRNG(), x0)
+    @test xx === Random.seed!(xx, xx)
+    @test TaskLocalRNG() === Random.seed!(TaskLocalRNG(), TaskLocalRNG())
+    @test xx == TaskLocalRNG()
+    @test xx == x1
+
+    # arrays
+    y2 = fork(xx, 2)
+    @test y2 isa Vector{Xoshiro} && size(y2) == (2,)
+    y34 = fork(xx, 3, 0x4)
+    @test y34 isa Matrix{Xoshiro} && size(y34) == (3, 4)
+    y123 = fork(xx, (1, 2, 3))
+    @test y123 isa Array{Xoshiro, 3} && size(y123) == (1, 2, 3)
+    y0 = fork(xx, ())
+    @test y0 isa Array{Xoshiro, 0} && size(y0) == ()
+    @test allunique([y2..., y34..., y123..., y0...])
+
+    # fork is currently unsupported for other RNGs
+    for rng = (RandomDevice(), MersenneTwister(0), TaskLocalRNG(), SeedHasher(0))
+        @test_throws MethodError fork(rng)
+        @test_throws MethodError fork(rng, (2, 3))
+    end
 end