julia version of mergeshuffle

mergeshuffle.jl

using Random: AbstractRNG, GLOBAL_RNG, MersenneTwister,shuffle!; import Future
using Base.Threads: nthreads, threadid
using BenchmarkTools

println("nthreads = $(nthreads())")

PAR_RNG = let m = MersenneTwister(1)
    [m; accumulate(Future.randjump, fill(big(10)^20, nthreads() - 1), init = m)]
end;


"""
    fisher_yates_shuffle!(v)
Fisher-Yates / Knuth shuffle algorithm
source code from https://www.rosettacode.org/wiki/fisher_yates__shuffle#Julia
"""
function fisher_yates_shuffle!(r::AbstractRNG, v::AbstractVector)
    @inbounds for i in length(v):-1:2
        j = rand(r, 1:i)
        v[i], v[j] = v[j], v[i]
    end
    return v
end
fisher_yates_shuffle!(v::AbstractVector) = fisher_yates_shuffle!(GLOBAL_RNG, v)

function merge!(r::AbstractRNG, v::AbstractVector, m::UInt, n::UInt)
    i::UInt = 1
    j::UInt = m + 1

    while true
        if rand(r, Bool)
            i == j && break
        else
            j == n && break
            @inbounds v[i], v[j] = v[j], v[i]
            j += 1
        end
        i += 1
    end
    @views shuffle!(r, v[i:n])
end
merge!(v::AbstractVector, m, n) = merge!(GLOBAL_RNG, v, m, n)

"""
    mergeshuffle!(v)
the paper: MergeShuffle: a very fast parallel random permutation algorithm
http://ceur-ws.org/Vol-2113/paper3.pdf
"""
function mergeshuffle!(r::Array{MersenneTwister,1}, v::AbstractVector; cutoff = 0x10000)
    n::UInt = length(v)
    # select q = 2^c such that n/q <= cutoff
    c::Int = 0
    while ((n >> c) > cutoff)
        c += 1
    end
    q::Int = 1 << c

    Threads.@threads for i in 0:(q - 1)
        j = (n * i) >> c
        k = (n * (i + 1)) >> c
        @views shuffle!(r[threadid()], v[j + 1:k])
    end
    p::UInt = 1
    while p < q
        Threads.@threads for i in 0:2p:(q - 1)
            j = (n * i) >> c
            k = (n * (i + p)) >> c
            l = (n * (i + 2p)) >> c
            @views merge!(r[threadid()], v[j + 1:l], k - j, l - j)
        end
        p += p
    end
    v
end
mergeshuffle!(v::AbstractVector) = mergeshuffle!([GLOBAL_RNG], v)



v = collect(1:10_000_000)
@time fisher_yates_shuffle!(v)
@time mergeshuffle!(PAR_RNG, v)
@time shuffle!(v)