Performance of threadsafe rand()

[mbauman]

I was surprised at how big a hit we suffer in looking up the THREAD_RNG — even when Julia is started with only 1 thread and without using @threads at all! Is there something we can do here?

julia> using Random, .Threads

julia> nthreads()
1

julia> function serialpi(n)
           inside = 0
           for i in 1:n
               x, y = rand(), rand()
               inside += (x^2 + y^2 <= 1)
           end
           return 4 * inside / n
       end
serialpi (generic function with 1 method)

julia> serialpi(1); @time serialpi(100_000_000);
  0.625890 seconds

julia> function serialpi2(n)
           inside = 0
           rng = Random.THREAD_RNGs[1]
           for i in 1:n
               x, y = rand(rng), rand(rng)
               inside += (x^2 + y^2 <= 1)
           end
           return 4 * inside / n
       end
serialpi2 (generic function with 1 method)

julia> serialpi2(1); @time serialpi2(100_000_000);
  0.257850 seconds

[KristofferC]

Ref #34216

[mbauman]

Ah, I knew we had this before — Kristoffer search is the best GitHub search.

Still, I feel like this is a slightly more constrained case; do we really need to be paying the cost of this indirection with nthreads()==1?

[tkf]

It'd be great if the compiler can hoist out the lookup for rng. It'd also help "context variables" as designed in #35833.

[yuyichao]

Ah, yeah, this is the reason I had before why thread local variable can't be replaced with task local and why caching of thread local variable can't be disabled. Unless there's a way to make RNG still work equally well when many more are created and destroied on the fly. Ref #35688

[mbauman]

This is now fixed due to #40546 (which removes the Random.THREAD_RNGs entirely). Equivalent test code, though, is:

julia> using Random, .Threads

julia> nthreads()
1

julia> function serialpi(n)
           inside = 0
           for i in 1:n
               x, y = rand(), rand()
               inside += (x^2 + y^2 <= 1)
           end
           return 4 * inside / n
       end
serialpi (generic function with 1 method)

julia> serialpi(1); @time serialpi(100_000_000);
  0.327586 seconds

julia> function serialpi2(n)
           inside = 0
           rng = Random.default_rng(1)
           for i in 1:n
               x, y = rand(rng), rand(rng)
               inside += (x^2 + y^2 <= 1)
           end
           return 4 * inside / n
       end
serialpi2 (generic function with 1 method)

julia> serialpi2(1); @time serialpi2(100_000_000);
  0.324798 seconds