A new RNG for Julia

gistfile1.jl

module XORShift

import Base.rand

type XORShiftStar1024 <: AbstractRNG
    p::Int
    state::Vector{Uint64}
end

XORShiftStar() = XORShiftStar1024(1, rand(Uint64, 16))

const globalXORShiftStar1024 = XORShiftStar()

function rand(rng::XORShiftStar1024)
    s = rng.state
    
    p = rng.p
    s0 = s[p]

    p = p %16 + 1
    s1 = s[p]

    s1 $= s1 << 31
    s1 $= s1 >> 11
    s0 $= s0 >> 30
    s0 $= s1
    s[p] = s0

    rng.p = p
    
    return s0 * 1181783497276652981
end

randXORShiftStar() = rand(globalXORShiftStar1024)*5.421010862427522e-20 # The number is 1/2.0^64

end #module

@stevengj Did you benchmark that? I would assume that the heterogeneous type results in type instability getfield (although you could always use unsafe_load if you wanted). Also I think the increments should be optimized out, since computing the address to load requires decrementing the argument to getfield/unsafe_load!.

I wonder whether this could be vectorized to compute two or four random numbers per iteration, but that would probably make getting a single number out of rand more complicated. I also don't think it can be made to work without llvmcall or maybe the SLP vectorizer.