add Random.jump(rng) API

We have long had methods for RNG "jumps ahead", i.e. advancing the state by
a given number of "steps", but no good API for that.

The only public API is `Future.randjump(r::MersenneTwister, steps::Integer)`,
and there are also functions for `Xoshiro` which are not public
(`Random.jump_128` and friends).

The following generic API is implemented here:
* `Random.jump(rng)` to jump by a reasonable default number of steps
* `Random.jump(rng; by::Real)` to jump by `by` steps
* `Random.jump!(rng; [by])` to equivalently jump in-place
* `Random.jump(rng, dims...; [by])` to create an array of jumped RNGs

In old julia versions, there also existed a method of `randjump` returning an
array, but the 1st element of this array was the passed argument;
the version here does not do this aliasing.

There are two kinds of integers one would wish to pass: dimensions for the array
version, and the number of steps.
Using jumps is relatively "niche", but needing to fidle with the number of steps
is even more niche. It's expected that in the vast majority of cases,
a good default is enough.

Some APIs in other languages have `jump` (e.g. 2^128 steps) and `long_jump`
(e.g. 2^192 steps), or `leap` in java, for more complicated cases;
for example each process gets a jumped RNG via
`long_jump`, and within each process, each thread gets a jumped RNG via `jump`.
But this is not very scalable if more kind of jumps are needed: should
`huge_jump` be introduced? For these rare cases where the default number of
steps is not sufficient, it seems better to let the programmer explicitly
specify the number of steps via an integer.

There is even a third kind of integers one might want to pass: in
`Random.jump_128(x::Xoshiro, i::Integer)`, `i` represents the number of times
a jump of size `2^128` is applied; this is because `Xoshiro` doesn't support
arbitrary number of steps; this is not supported in the proposed API, because
1) it's trivial for the user to implement herself, and 2) in probably most use
cases, using the array version will be a valid alternative, and more efficient
because previous computations are not wasted
(like in `[Random.jump_128(x, i) for i=1:num_tasks]` vs
`Random.jump(x, num_tasks)`).

Another argument in favor of this API is that it mirrors the
proposed `Random.fork(rng, dims...)` function from #58193.

Author: rfourquet

NEWS.md

@@ -47,6 +47,10 @@ Standard library changes
 
 #### Profile
 
+#### Random
+
+* New `Random.jump` function to advance the state ("jump ahead") of `Xoshiro` or `MersenneTwister` RNGs ([#58353]).
+
 #### REPL
 
 #### Test

stdlib/Future/src/Future.jl

@@ -35,11 +35,10 @@ Create an initialized `MersenneTwister` object, whose state is moved forward
 One such step corresponds to the generation of two `Float64` numbers.
 For each different value of `steps`, a large polynomial has to be generated internally.
 One is already pre-computed for `steps=big(10)^20`.
+
+!!! compat "Julia 1.13"
+    As of Julia 1.13, this functionality is now implemented by `Random.jump`.
 """
-function randjump(r::MersenneTwister, steps::Integer)
-    j = Random._randjump(r, Random.DSFMT.calc_jump(steps))
-    j.adv_jump += 2*big(steps) # convert to BigInt to prevent overflow
-    j
-end
+randjump(r::MersenneTwister, steps::Integer) = Random.jump(r; by=steps)
 
 end # module Future

stdlib/Random/docs/src/index.md

@@ -83,6 +83,8 @@ Random.TaskLocalRNG
 Random.Xoshiro
 Random.MersenneTwister
 Random.RandomDevice
+Random.jump
+Random.jump!
 ```
 
 ## [Hooking into the `Random` API](@id rand-api-hook)

stdlib/Random/src/MersenneTwister.jl

@@ -51,6 +51,14 @@ The `seed` may be an integer, a string, or a vector of `UInt32` integers.
 If no seed is provided, a randomly generated one is created (using entropy from the system).
 See the [`seed!`](@ref) function for reseeding an already existing `MersenneTwister` object.
 
+`MersenneTwister` supports "jumping ahead" via [`Random.jump`](@ref), with an
+arbitrary number of steps, where one step represents consuming two `Float64` values.
+Given `MersenneTwister`'s huge period of `2^19937-1`, there is a lot of flexibility
+for arranging -- via different numbers of jump steps -- computations with multiple
+levels of parallelism requiring non-overlapping random subsequences.
+Prior to Julia 1.13, the same functionality was available from the `Future` stdlib
+with the `Future.randjump` function.
+
 !!! compat "Julia 1.11"
     Passing a negative integer seed requires at least Julia 1.11.
 
@@ -554,7 +562,7 @@ function rand!(r::MersenneTwister, A1::Array{Bool}, sp::SamplerType{Bool})
 end
 
 
-### randjump
+### jump
 
 # Old randjump methods are deprecated, the scalar version is in the Future module.
 
@@ -568,17 +576,16 @@ function _randjump(r::MersenneTwister, jumppoly::DSFMT.GF2X)
     s
 end
 
-# NON-PUBLIC
-function jump(r::MersenneTwister, steps::Integer)
-    iseven(steps) || throw(DomainError(steps, "steps must be even"))
-    # steps >= 0 checked in calc_jump (`steps >> 1 < 0` if `steps < 0`)
-    j = _randjump(r, Random.DSFMT.calc_jump(steps >> 1))
-    j.adv_jump += steps
-    j
+function jump(rng::MersenneTwister; by::Real=NaN)
+    isnan(by) && (by = 2.0^128)
+    steps = BigInt(by)
+    # steps >= 0 checked in calc_jump
+    jumped = _randjump(rng, DSFMT.calc_jump(steps))
+    jumped.adv_jump += steps
+    jumped
 end
 
-# NON-PUBLIC
-jump!(r::MersenneTwister, steps::Integer) = copy!(r, jump(r, steps))
+jump!(rng::MersenneTwister; by::Real=NaN) = copy!(rng, jump(rng; by))
 
 
 ### constructors matching show (EXPERIMENTAL)
@@ -642,7 +649,7 @@ function advance!(r::MersenneTwister, adv_jump, adv, adv_vals, idxF, adv_ints, i
     ms = dsfmt_get_min_array_size() % Int
     work = sizehint!(Vector{Float64}(), 2ms)
 
-    adv_jump != 0 && jump!(r, adv_jump)
+    adv_jump != 0 && jump!(r; by=adv_jump)
     advF = (adv_vals, idxF) != (0, 0)
     advI = (adv_ints, idxI) != (0, 0)
 

stdlib/Random/src/RNGs.jl

@@ -306,3 +306,92 @@ function hash_seed(str::AbstractString, ctx::SHA_CTX)
     end
     SHA.update!(ctx, (0x05,))
 end
+
+
+## jump
+
+"""
+    Random.jump(rng::AbstractRNG, [dims...]; [by::Real])
+
+Create a new generator of the same type as `rng`, whose state is advanced forward
+("jump ahead") by `by` "steps".
+This is equivalent to `Random.jump!(copy(rng); by)`.
+
+When `dims` is specified (as integers or as a tuple), create an array of such
+generators `x_1, x_2, ..., x_n`, where `x_1 = jump(rng; by)`, and
+`x_i = jump(x_{i-1}; by)`.
+
+This can be used to generate non-overlapping subsequences for parallel computations.
+The maximum number of such subsequences is roughly equal to the period
+of `rng` divided by `by`.
+
+The default value of `by` should generally be large enough such that the generated
+subsequences can accomodate any computation;
+it's currently `2.0^128` for `MersenneTwister` and `Xoshiro`, but can change in the
+future.
+
+!!! compat "Julia 1.13"
+    This function was introduced in Julia 1.13.
+
+# Examples
+```julia-repl
+julia> Random.jump(Xoshiro(0))
+Xoshiro(0xe223c163bccd575d, 0xef13c6b1ba8b980b, 0xeb86b37bd43e78ab, 0x2b93c9c43b0815fb, 0x22a21880af5dc689)
+
+julia> m = MersenneTwister(123)
+MersenneTwister(123)
+
+julia> Random.jump(m, 3; by=1000)
+3-element Vector{MersenneTwister}:
+ MersenneTwister(123, (1000, 0))
+ MersenneTwister(123, (2000, 0))
+ MersenneTwister(123, (3000, 0))
+```
+
+See also [`Random.jump!`](@ref).
+"""
+function jump end
+
+"""
+    Random.jump!(rng::AbstractRNG; [by::Real]) -> rng
+
+Advance the state of `rng` forward ("jump ahead") by `by` "steps".
+What a step represents is specific to each specific type.
+For example with `Xoshiro`, advancing by one step means consuming 8 bytes from `rng`.
+
+Note that some approximation in the number of steps is permitted in the implementation.
+For example, if the RNG caches some random numbers, it's allowed to discard the cache
+and jump ahead by `by` steps without accounting for the specific number of discarded
+values. In general, `by` should be specified large enough, with a huge margin, such
+that such details do not matter.
+
+`by` should be an integer, but can be expressed via non-`Integer` types for
+convenience, e.g. `by = 2.0^128`. While some RNGs like `MersenneTwister` support any
+positive integer, others like `Xoshiro` support only specific values.
+
+!!! compat "Julia 1.13"
+    This function was introduced in Julia 1.13.
+
+# Examples
+```julia-repl
+julia> m = MersenneTwister(123)
+MersenneTwister(123)
+
+julia> Random.jump!(m; by=1000)
+MersenneTwister(123, (1000, 0))
+```
+See also [`Random.jump`](@ref).
+"""
+function jump! end
+
+jump(rng::AbstractRNG; by::Real=NaN) = jump!(copy(rng); by)
+jump(rng::AbstractRNG, d1::Integer, dims::Integer...; by::Real=NaN) =
+    jump(rng, Dims((d1, dims...)); by)
+
+function jump(rng::AbstractRNG, dims::Dims; by::Real=NaN)
+    js = Array{typeof(rng)}(undef, dims)
+    for ii in eachindex(js)
+        rng = js[ii] = jump(rng; by)
+    end
+    js
+end

stdlib/Random/src/Random.jl

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