task splitting: change additive accumulation to multiplicative

[StefanKarpinski]

Issue raised here. Given this definition, look at output values:

using .Threads

macro fetch(ex) :(fetch(@spawn($(esc(ex))))) end

function taskCorrelatedXoshiros(+ = +)
    r11, r10 = @fetch (@fetch(rand(UInt64)), rand(UInt64))
    r01, r00 = (@fetch(rand(UInt64)), rand(UInt64))
    r01 + r10 - (r00 + r11)
end

Before:

julia> sort!(unique(taskCorrelatedXoshiros() for _ = 1:1000))
9-element Vector{UInt64}:
 0x0000000000000000
 0x0000000000000001
 0x00000000007fffff
 0x0000000000800000
 0x0000000000800001
 0xffffffffff7fffff
 0xffffffffff800000
 0xffffffffff800001
 0xffffffffffffffff

After:

julia> sort!(unique(taskCorrelatedXoshiros() for _ = 1:1000))
1000-element Vector{UInt64}:
 0x000420647a085960
 0x0038c5b889b585c6
 0x007b29fae8107ff7
 0x00e73b9e883ac1c8
                  ⋮
 0xfe68be9c0dde1e88
 0xfeca042354218c35
 0xfeeb8203e470c96b
 0xfff5dbb8771307b9

julia> sort!(unique(taskCorrelatedXoshiros(*) for _ = 1:1000))
1000-element Vector{UInt64}:
 0x00126f951c1e56dc
 0x0025a82477ffac08
 0x002dd82c9986457a
 0x00a713c4d56a3dbc
                  ⋮
 0xfe2e40a5b345095e
 0xfe77b90881967436
 0xfea2559be63f1701
 0xff88b5b28cefac5f

[vtjnash]

Where do these new a[1:4] values come from?

[Seelengrab]

I guess this part of the explanation/doc above jl_rng_split should be adjusted too:

It also means that when accumulating the dot product incrementally, as
described in SplitMix, we don't need to multiply weights by anything, we simply
add the random weight for the current task tree depth to the parent's dot
product to derive the child's dot product.

[StefanKarpinski]

Where do these new a[1:4] values come from?

As the comment says, they are "random additive offsets".

[StefanKarpinski]

I guess this part of the explanation/doc above jl_rng_split should be adjusted too:

Yes, I'm not adjusting and comment text until we settle on something.

[oscardssmith]

Does this change notably affect the task spawn time? (I'm assuming not, because add, multiply and byteswap are all really fast, but good to make sure). Otherwise, this looks like a good change.

[StefanKarpinski]

I get identical timings between master and this rebased onto master—about 65ns. Would be good if someone could benchmark on different systems. I'm on an M2 system, so x86 would be helpful.

[carstenbauer]

Just ran a quick and dirty benchmark on one of our compute nodes (x86 AMD Zen 3).

# master (a1db8daafd*)
julia> @btime @spawn(nothing);
  65.741 ns (4 allocations: 480 bytes)

# this PR rebased onto master
julia> @btime @spawn(nothing);
  66.098 ns (4 allocations: 480 bytes)

[StefanKarpinski]

Ok, this is ready to go from my perspective.

[KristofferC]

I want to see what package's tests "broke" from the RNG stream changing here.

@nanosoldier runtests(; vs = "@c80a9647166d59782eff6be8c75a5052e17483d7")

[StefanKarpinski]

Hopefully not too disruptive since it only changes the RNG stream in child tasks, not the seeded task; I suspect that it's fairly rare for packages to do multithreaded deterministic RNG testing.