How can I use randjump instead of generating a certain number of random numbers?

user_7553531 0
Viewed 80

I am trying to "skip forward" a few realizations by using the function Future.randjump(), but it doesn't seem to behave as I expect it to. The following code gives me the desired result, where jumping forward 1 steps gives the same result as if I had called rand(rng) twice, i.e. the two println display the same number:

using Random, Future
rng = MersenneTwister(123);

new_rng = Future.randjump(rng, 1)
rand(rng)
rand(rng)
println(rand(rng))
println(rand(new_rng))

However, if I add one extra call to rand(rng) before the call to randjump(), the two printed numbers are completely different:

using Random, Future
rng = MersenneTwister(123);

rand(rng)                                # Added line
new_rng = Future.randjump(rng, 1)
rand(rng)
rand(rng)
println(rand(rng))
println(rand(new_rng))

I expected that the two calls to println() would display the same thing even in the second case, how come they don't? Is there a way I can use randjump() in the second case to get the same realizations as if I had called rand(rng) several times? Thank you in advance.

julia random
1 Answers

One unit of randjump corresponds to generation of two floating point numbers.

Consider this example

julia> rng = MersenneTwister(123);

julia> rng2 = Future.randjump(rng, 1);

julia> rand(rng, 4)
4-element Vector{Float64}:
 0.7684476751965699
 0.940515000715187
 0.6739586945680673
 0.3954531123351086

julia> rand(rng2,2)
2-element Vector{Float64}:
 0.6739586945680673
 0.3954531123351086

Note that in the second call (that is rand(rng2,2)) the both numbers are identical to the two last numbers in the first call (taht is rand(rng,2)).

Another issue is that different distributions might "consume" Float64 numbers from the stream at a different speed - so you need to check with a particular distribution how fast it consumes floats for the stream (some might also use buffering etc...).

Looking at the source code of randn (@edit randn()) it consumes one float and hence you get the same results for those two calls:

julia> randn(MersenneTwister(123),6)[3:end]
4-element Vector{Float64}:
  1.142650902867199
  0.45941562040708034
 -0.396679079295223
 -0.6647125451916877

julia> randn(Future.randjump(MersenneTwister(123),1),4)
4-element Vector{Float64}:
  1.142650902867199
  0.45941562040708034
 -0.396679079295223
 -0.6647125451916877

EDIT

Regarding your comment the size of Mersenne Twister state is 19937 bits and half-unit jumps are not supported. Running rand is mutating this state but not half-the way - so you end up with different bits. Note that an RNG is a sequence of states and the actual values are calculated from that state.

The correct pattern to synchronize random numbers in your computations is the following:

master_rng = MersenneTwister(123);

rng1 = Future.randjump(master_rng, big(10)^20)
# do whatever you want

rng2 = Future.randjump(master_rng, 2*big(10)^20)
# do whatever you want

rng3 = Future.randjump(master_rng, 3*big(10)^20)
# do whatever you want

With this pattern you can correctly maintains synchronization between random number streams and have full control whether the should overlap or not.

user_99577100
Related