Implement cong in Julia

[vchuravy]

Moving cong to Julia safes a tiny bit in overhead,
but accessing the PTLS rngseed for the scheduler is a bit... ugly.

[gbaraldi]

Can we do something that doesn't need the PTLS?

[vchuravy]

I think the challenge is that we want randomness that is independent of the task state, since we don't have "scheduler state" outside of the ptls that's the only option.

I am more annoyed how slow this is due to the mod operation.

[oscardssmith]

Can we guarantee that the size of the heaps are powers of 2? If so, you can just use an &.

[vchuravy]

It's the number of heaps and not the size of them, and that is equal to the number of threads. Which often is a power of two, but not guaranteed

[oscardssmith]

What about using Base.MultiplicativeInverses to precompute the multiplicativeinverse of the number of threads?

[IanButterworth]

Is this a valid benchmark?

Update: ** not a fair test because the base commit is different **

this PR

% ./julia -t6 -q
julia> function fib(n::Int)
                  if n < 2
                      return n
                  end
                  t = Threads.@spawn fib(n - 2)
                  return fib(n - 1) + fetch(t)
              end
fib (generic function with 1 method)

julia> using BenchmarkTools

julia> @btime fib(30)
  2.098 s (6929222 allocations: 1.13 GiB)
832040

julia> @btime fib(30)
  2.029 s (6896970 allocations: 1.05 GiB)
832040

master b88f64f (via juliaup)

% julia +nightly -t6
...
julia> @btime fib(30)
  748.326 ms (6772487 allocations: 734.01 MiB)
832040

julia> @btime fib(30)
  700.210 ms (6743099 allocations: 659.97 MiB)
832040

[vchuravy]

Probably. I haven't checked the type-stability of this code after the recent rewrite. Since the allocation increased I assume we have one.

@IanButterworth do you want to take this over?

[vchuravy]

@vtjnash any ideas on how to do this nicer?

[IanButterworth]

I made it so we only unsafe_convert plts once. Is that nice enough?

[IanButterworth]

Doesn't appear to be an allocation within cong. And it's faster.
This PR b5d065a

julia> @btime Base.Partr.cong(UInt32(100))
  2.333 ns (0 allocations: 0 bytes)

master

julia> @btime Base.Partr.cong(UInt32(100))
  5.333 ns (0 allocations: 0 bytes)
0x00000036

However I just locally updated this branch from master and it regressed to

julia> @btime Base.Partr.cong(UInt32(100))
  6.958 ns (0 allocations: 0 bytes)
0x00000061

[gbaraldi]

So looking into this together with Ian we discovered that our rng is just very bad. I also discovered that get_ptls is not getting optimized as it should because the declaration is behind an ifdef so the ccall transform doesn't do it.
Finally we also discovered that rand(1:range) performs quite on par with the current cong and better than the corrected one.
The corrected one being
Following the origin in https://opensource.apple.com/source/Libc/Libc-1439.141.1/gen/FreeBSD/arc4random.3.auto.html
instead of the wrong implementation in https://www.pcg-random.org/posts/bounded-rands.html.
Basically the current thing only uses n bits where n is 2^n for minimal value that's largest than the max. Which isn't a lot of bits leading to basically cycling over the numbers.

function _cong(max::UInt64, seed::UInt64)
    if max < 2
        return UInt64(0), seed
    end
    mask = typemax(UInt64)
    zeros = leading_zeros(max)
    mask >>= zeros
    bits = 8 * sizeof(UInt64) - zeros
    while true
        value = UInt64(69069) * seed + UInt64(362437)
        seed = value
        x = seed & mask
        if (x < max)
            return x, seed
        end
        bits_left = zeros
        while bits_left >= bits
            value >>= bits
            x = value & mask
            if x < max
                return x, seed
            end
            bits_left -= bits
        end
    end
    return x, seed
end

[IanButterworth]

To illustrate the periodicity of cong on master.

julia> for i in 2:20
           @info i
           x = [Base.Partr.cong(UInt32(i)) for _ in 1:100];
           display(lineplot(x, width=100, height = 8))
       end

2-6	7-11	12-16