Is `rand()` threadsafe?

While it seems to work and the output looks sensible, I have this lingering suspicion that rand() is not threadsafe and that using it in a Thread.@threads loop is volatile.

Browsing through base/random.jl and base/dSFMT.jl I saw no obvious sign of thread protection, which makes me wonder why it seems to work so well in the first place. Is dSFMT itself threadsafe maybe? I did not find any info on that.

So my question is if it is safe to call rand() in a Threads.@threads loop, and if not what the encourage way of doing so is.

Something like this?

julia> const trand_mutex = Threads.Mutex()
Base.Threads.Mutex(0, Ptr{Void} @0x00007ff8949c3fc0)

julia> function trand(args...)
           lock(trand_mutex)
           result = rand(args...)
           unlock(trand_mutex)
           result
       end
trand (generic function with 1 method)

Is rand() threadsafe?

No. It won’t crash since there’s nothing really for it to crash, but it’ll not return the correct result and perform poorly.

Create a rng per thread. What I’ve been doing is here

Definitely not.

1 Like

Thanks for the quick reply!

I know very little about random number generators, but I was under the impression that it is generally discouraged to use multiple RNG in parallel.

from A Primer on Repeatable Random Numbers | Unity Blog

No it just say that you need to be careful about initializing it. For my use case, I don’t care.

Understood, thanks for the help and example!

Won’t this make every thread have identical random numbers?

@noinline function init_thread_rng()
    # Allocate the random number generator on the thread's own heap lazily
    # instead of the master thread heap to minimize memory conflict.
    ThreadRNG[Threads.threadid()] = MersenneTwister(0)
end

I assume for separate random numbers per thread, it’d be fine to initialize with MersenneTwister(Threads.threadid())?

Yes.

My understanding of seeding is that method (of using consecutive small integers) isn’t necessarily any random than MersenneTwister(0), but that you can jump one RNG ahead to effectively split it into two (or N) independent streams (https://github.com/JuliaLang/julia/pull/12498). I could be wrong though.

1 Like

Hi all

random.jl has a test which shows the usage of randjump
https://github.com/JuliaLang/julia/blob/c1f107fa99939ede9f963fef4044097215c77977/test/random.jl#L468

I understand that the default jump in julia is 10^20.
So if I need 36 random generators for 36 processes and on each thread I anticipate to generate around 10^13 random numbers. Am I safe by using the code below?
Will each MersennTwister be 10^20 ‘apart’ from the next one? (the documentation is not too easy to read)?

seed=2001
mta = MersenneTwister(seed)
mts = randjump(mta, 36)
# -> use mts[x] on process x (for x in 1:36), rand(mts[x]) respectively 

If anyone could provide a jump polynomial for higher steps that would be nice. I was not able to use the code from the site below (as my C skills are nil)
http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/JUMP/dsfmt-jump.html

Thank you

Am I safe by using the code below?

yes

Will each MersennTwister be 10^20 ‘apart’ from the next one?

yes

If anyone could provide a jump polynomial for higher steps that would be nice

You can use the code at MersenneTwister: implement generation of jump-ahead polynomials by rfourquet · Pull Request #16906 · JuliaLang/julia · GitHub (which bascically is a port from the C code you mention) and why not give support to this PR to help get it merged!

Thanks for the quick answer.

I now realize that I may run into issues when generating that man parallel ‘parts’ of a single mersenneTwister (considering its period of 10^61)

However it seems that there is apackage which solves my issues and has ‘better’ rngs:

http://sunoru.github.io/RandomNumbers.jl/latest/

Are you sure?
Not 2^19937?

Oh. My bad.
This is embarrassing, I must have confused 10^6001 (which should be close to 2^19937) with 10^61.

I guess this can only happen to someone with no intuition for power series.
You are correct.
Seems like I have plenty of ‘room’ available with 36*20=720 << 6001

That is to say: It is my understanding that the second element of the array would effectively start at the ‘random element’ 1+10^20 which the first MersenneTwister of the array produces.

Just in case anyone stumbles on this old thread, rand() and friends are now thread-safe (as of v1.3). See here for more info.

5 Likes