# Random numbers and threads

**URL:** <https://discourse.julialang.org/t/random-numbers-and-threads/77364>\
**Category:** General Usage\
**Tags:** question, multithreading, random\
**Created:** [March 3, 2022, 3:31pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364 "2022-03-03T15:31:24Z")\
**Posts on this page:** 11\
**Page:** 2

<div class="post-metadata">

**Author:** ![tkf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkf/32/17635_2.png) [@tkf](https://discourse.julialang.org/u/tkf)\
**Post date:** [March 4, 2022, 6:45am UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/21 "2022-03-04T06:45:45Z")

</div>

> [@jw3126](#):
>
> If I spawn things in my main task, its rng state changes.

This is almost correct, but it actually happens during `Task` creation (as @fredrikekre has already mentioned)

```julia
julia> Task(nothing)
       Random.seed!(0)
       rand()
0.4056994708920292

julia> Random.seed!(0)
       Task(nothing)
       rand()
0.6616126907308237

```

> [@jw3126](#):
>
> As of julia 1.7

Note that, in Julia 1.7, `@threads`’s scheduling is context-dependent and thus the task tree is not quite predictable. If you use it `@threads` in, e.g., at the “root task,” you get multiple tasks (on multiple worker threads):

```julia
julia> Threads.@threads for j in 1:Threads.nthreads()
           @show j => current_task()
       end
j => current_task() = 2 => Task (runnable) @0x00007fb46da0dcd0
j => current_task() = 4 => Task (runnable) @0x00007fb46da0dfb0
j => current_task() = 3 => Task (runnable) @0x00007fb46da0de40
j => current_task() = 1 => Task (runnable) @0x00007fb46da0db60

```

However, inside of another `@threads`, it uses single task (on single worker thread):

```julia
julia> Threads.@threads for i in 1:Threads.nthreads()
           if i == 1
               Threads.@threads for j in 1:Threads.nthreads()
                   @show j => current_task()
               end
           end
       end
j => current_task() = 1 => Task (runnable) @0x00007fb46e331880
j => current_task() = 2 => Task (runnable) @0x00007fb46e331880
j => current_task() = 3 => Task (runnable) @0x00007fb46e331880
j => current_task() = 4 => Task (runnable) @0x00007fb46e331880

```

Thus, it implies that `Random.seed!(0); f()` may not have the same result, depending on the context. This is fixed in Julia 1.8 thanks to PRs such as

> <https://github.com/JuliaLang/julia/pull/43919>
>
> Ref https://github.com/JuliaLang/julia/pull/35646 - the idea here is largely pro…posed there. This is perhaps mainly a question of whether now is the right time to do this?
> 
> \_\_\_
> 
> Enables a more composable variant of \`Threads.@threads\` via a new \`:dynamic\` schedule strategy.
> 
> An example where \`busywait\` is a non-yielding timed loop that waits for a number of seconds.
> Note that \`:static\` is the current default where no schedule is specified:
> 
> \`\`\`julia
> julia\> @time begin
> Threads.@spawn busywait(5)
> Threads.@threads :static for i in 1:Threads.nthreads()
> busywait(1)
> end
> end
> 6.003001 seconds (16.33 k allocations: 899.255 KiB, 0.25% compilation time)
> 
> julia\> @time begin
> Threads.@spawn busywait(5)
> Threads.@threads :dynamic for i in 1:Threads.nthreads()
> busywait(1)
> end
> end
> 2.012056 seconds (16.05 k allocations: 883.919 KiB, 0.66% compilation time)
> \`\`\`
> 
> Unlike \`:static\` the \`:dynamic\` strategy can be nested.
> 
> \`\`\`julia
> julia\> @time begin
> Threads.@spawn busywait(5)
> Threads.@threads :dynamic for i in 1:Threads.nthreads()
> Threads.@threads :dynamic for j in 1:Threads.nthreads()
> busywait(1)
> end
> end
> end
> 17.023767 seconds (61.46 k allocations: 3.402 MiB, 0.14% compilation time)
> 
> julia\> Threads.nthreads()
> 16
> \`\`\`
> 
> @vtjnash: A thought was that the threaded region didn't need to be marked with this strategy, is that correct?
> 
> 
> Co-authors: @tkf @jpsamaroo

As for manual task creation that @lmiq is mentioning, I think the most important benefit is that it lets you have multiple streams of RNG and also with arbitrary algorithms.

---

<div class="post-metadata">

**Author:** ![jw3126](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jw3126/32/3086_2.png) [@jw3126](https://discourse.julialang.org/u/jw3126)\
**Post date:** [March 4, 2022, 8:15am UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/22 "2022-03-04T08:15:32Z")

</div>

Thanks a lot, @tkf @lmiq @fredrikekre @stevengj I now understand the global rng and how I could use the “task pool pattern” mentioned by @lmiq . I have one more question. I used the pattern in the opening post

```julia
rngs = [make_rng(tid) for tid in 1:Threads.nthreads()]

```

quite a bit in old code that I don’t want to rewrite.

- Is this pattern safe in Julia 1.7?
- Is there a super simple way to make it safe? Say by annotating each loop like `@threads :static for ...` or something?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 4, 2022, 1:39pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/23 "2022-03-04T13:39:59Z")

</div>

> [@jw3126](#):
>
> Is this pattern safe in Julia 1.7?

No, as I [commented above](https://discourse.julialang.org/t/random-numbers-and-threads/77364/8).

You could use task-local storage, e.g. something like:

```julia
rng = get!(task_local_storage(), _RNGS_KEY) do
    MersenneTwister()
end::MersenneTwister

```

where `const _RNGS_KEY = gensym(:rngs)`

But since I think new tasks are spawned every time you do `@threads`, it’s not clear what the advantage of this would be over just `rng = MersenneTwister()` in the type of usage you seem like you’re talking about.

---

<div class="post-metadata">

**Author:** ![tkf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkf/32/17635_2.png) [@tkf](https://discourse.julialang.org/u/tkf)\
**Post date:** [March 4, 2022, 2:34pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/24 "2022-03-04T14:34:09Z")

</div>

> [@jw3126](#):
>
> annotating each loop like `@threads :static for ...`

I actually recommend this, if you just want a very easy and conservative approach. There’s no change in Julia 1.7 about this but it could have been invalid already in 1.5 and later if you happen to start using old code in nested tasks. Furthermore, there are changes in the upcoming 1.8 which makes old programs invalid in more cases.

For more information, we’ve updated the docstring to explain the semantics for Julia 1.8 [Multi-Threading · The Julia Language](https://docs.julialang.org/en/v1.9-dev/base/multi-threading/#Base.Threads.@threads)

But the short answer is that this is because the following code may throw in 1.8 and later

```julia
Threads.@threads for _ in 1:Threads.nthreads()
    i = Threads.threadid()
    yield()
    j = Threads.threadid()
    @assert i == j
end

```

> [@jw3126](#):
>
> ```julia
> rngs = [make_rng(tid) for tid in 1:Threads.nthreads()]
> 
> ```

The problem is rather how this is used. If you access it using “task id” you issued yourself, using a pattern like [Random numbers and threads - #20 by lmiq](https://discourse.julialang.org/t/random-numbers-and-threads/77364/20), there’s no problem. But if you use `threadid()`, there is a problem in 1.8 or later as I mentioned above. If you use `rngs` as `rand(rngs[threadid()])` or something similar, it may be fine (but it depends on the RNG type implementation). But, if you cache `rngs` access in a variable, it’d be problematic:

```julia
Threads.@threads for _ in 1:Threads.nthreads()
    rng = rngs[Threads.threadid()]
    f() # some function that may contain a yield point (e.g., printing, locking, ...)
    rand(rng) # not OK because we may have `rng !== rngs[Threads.threadid()]`
end

```

The rule of thumb for good high-level multi-tasking code in Julia is “don’t talk about threads” although this attitude is rather too idealistic ATM.

---

<div class="post-metadata">

**Author:** ![Red-Portal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/red-portal/32/9102_2.png) [@Red-Portal](https://discourse.julialang.org/u/Red-Portal)\
**Post date:** [March 6, 2022, 10:24pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/25 "2022-03-06T22:24:30Z")

</div>

Hi,

I would like to mention that you shouldn’t combine the random number sequence of Mersenne Twister RNGs initialized with different seeds. The theoretical guarantees of Mersenne Twisters don’t work in those settings. Instead, it is necessary to use RNGs specialized for parallel RNG.

Some typical choices for parallel RNG are the [PCG](https://www.pcg-random.org/index.html) and [Random123](https://ieeexplore.ieee.org/document/6114424) families. These RNGs work by having another layer of “seeds” known as keys. You start from a single seed, and generate multiple RNG sequences by passing a key for each parallel sequence.

While both are not “cryptographically secure”, they both generate good quality random numbers with much less computation compared to Mersenne Twisters. They are also quite simple to use. Personally, I prefer Random123, and there is a nice [implementation](https://github.com/JuliaRandom/Random123.jl) in Julia.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [March 6, 2022, 10:44pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/26 "2022-03-06T22:44:58Z")

</div>

I beg your pardon. Are you saying that Julia Base’s random number generators aren’t thread safe? I thought this had all been worked out.

---

<div class="post-metadata">

**Author:** ![Red-Portal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/red-portal/32/9102_2.png) [@Red-Portal](https://discourse.julialang.org/u/Red-Portal)\
**Post date:** [March 6, 2022, 10:59pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/27 "2022-03-06T22:59:38Z")

</div>

Hi,

No, what I’m saying is that you shouldn’t use Mersenne Twisters in parallel even if they are thread safe. It’s “mathematically” unsafe. This is a very common mistake that is even stated in [Wikipedia](https://en.wikipedia.org/wiki/Mersenne_Twister#Disadvantages):

> Multiple Mersenne Twister instances that differ only in seed value (but not other parameters) are not generally appropriate for Monte-Carlo simulations that require independent random number generators, though there exists a method for choosing multiple sets of parameter values.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [March 6, 2022, 11:15pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/28 "2022-03-06T23:15:24Z")

</div>

Is this peculiar to Mersenne Twisters? I understand that they are no longer the default rng in Julia, so perhaps that is now less of an issue?

---

<div class="post-metadata">

**Author:** ![Red-Portal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/red-portal/32/9102_2.png) [@Red-Portal](https://discourse.julialang.org/u/Red-Portal)\
**Post date:** [March 6, 2022, 11:24pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/29 "2022-03-06T23:24:54Z")

</div>

No, it’s a problem with any RNG that is not specifically designed to be used in parallel. To do proper parallel RNG generation, you _have_ to use RNGs specialized for parallel generation (like Random123) in the correct way.

---

<div class="post-metadata">

**Author:** ![Sukera](https://avatars.discourse-cdn.com/v4/letter/s/ce7236/32.png) [@Sukera](https://discourse.julialang.org/u/Sukera)\
**Post date:** [March 29, 2022, 12:21pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/30 "2022-03-29T12:21:05Z")

</div>

Are the arguments surrounding Mersenne Twister even applicable to Xoshiro? According to the authors, [they do provide jump polynomials](https://prng.di.unimi.it) for this exact use case:

> All generators, being based on linear recurrences, provide _jump functions_ that make it possible to simulate any number of calls to the next-state function in constant time, once a suitable _jump polynomial_ has been computed. We provide ready-made jump functions for a number of calls equal to the square root of the period, to make it easy generating non-overlapping sequences for parallel computations, and equal to the cube of the fourth root of the period, to make it possible to generate independent sequences on different parallel processors.

Though admittedly, I don’t think they are exposed in julia. Maybe they should be added to the `Random` stdlib.

---

<div class="post-metadata">

**Author:** ![Palli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/palli/32/3380_2.png) [@Palli](https://discourse.julialang.org/u/Palli)\
**Post date:** [March 29, 2022, 2:31pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/31 "2022-03-29T14:31:53Z")

</div>

> [@lmiq](#):
>
> Xoshiro is slower in that minimal test, effectively

That comes as a surprise (something that can be optimized, or at least faster in all real-word code?). Xoshiro supposed to be fast, it’s short code, with much smaller state than MersenneTwister.

> [@stevengj](#):
>
> the “global” rng is actually a task-local state.

Right, as of 1.7.

> [@stevengj](#):
>
> > [@lmiq](#):
> >
> > `rng = Random.MersenneTwister()`
> 
> This is not even the same rng as the default (Xoshiro).

If you mention Xoshiro by name in your code, it’s no longer compatible with Julia 1.6 LTS. If you just use a global thread, I suppose you can have the same code working in both 1.7/1.8 and 1.6 LTS, but might there be a way for 1.7 code using threads to keep some (partial) compatibility with 1.6? In both cases the stream of numbers will be different.

[If you want to target 1.6 and later, then you can use MersenneTwister, that’s I believe why it was kept around, for compatibility, for those who mentioned it my name (thus will get same stream).]

[Previous page](https://discourse.julialang.org/t/random-numbers-and-threads/77364.md?page=1)