CODEWITHSUNDEEP
randomintegerjulia
epx
epx

Reputation: 809

Efficient way of generating random integers within a range in Julia

I'm doing MC simulations and I need to generate random integers within a range between 1 and a variable upper limit n_mol

The specific Julia function for doing this is rand(1:n_mol) where n_mol is an integer that changes with every MC iteration. The problem is that doing it this is slow... (possibly an issue to open for Julia developers). So, instead of using that particular function call, I thought about generating a random float in [0,1) multiply it by n_mol and then get the integer part of the result: int(rand()*n_mol) the problem now is that int() rounds up so I could end up with numbers between 0 and n_mol, and I can't get 0... so the solution I'm using for the moment is using ifloor and add a 1, ifloor(rand()*n_mol)+1, which considerably faster that the first, but slower than the second.

function t1(N,n_mol)
    for i = 1:N
        rand(1:n_mol)
    end
end

function t2(N,n_mol)
    for i = 1:N
        int(rand()*n_mol)
    end
end

function t3(N,n_mol)
    for i = 1:N
        ifloor(rand()*n_mol)+1
    end
end


@time t1(1e8,123456789)
@time t2(1e8,123456789)
@time t3(1e8,123456789)


elapsed time: 3.256220849 seconds (176 bytes allocated)
elapsed time: 0.482307467 seconds (176 bytes allocated)
elapsed time: 0.975422095 seconds (176 bytes allocated)

So, is there any way of doing this faster with speeds near the second test? It's important because the MC simulation goes for more than 1e10 iterations. The result has to be an integer because it will be used as an index of an array.

Upvotes: 4

Views: 2639

Answers (2)

Moritz Schauer
Moritz Schauer

Reputation: 1417

The rand(r::Range) code is quite fast, given the following two considerations. First, julia calls a 52 bit rng twice to obtain random integers and a 52 bit rng once to obtain random floats, that gives with some book keeping a factor 2.5. A second thing is that 

(rand(Uint) % k)

is only evenly distributed between 0 to k-1, if k is a power of 2. This is taken care of with rejection sampling, this explains more or less the remaining additional cost.

If speed is extremely important you can use a simpler random number generator as Julia and ignore those issues. For example with a linear congruential generator without rejection sampling 

function lcg(old) 
    a = unsigned(2862933555777941757)
    b = unsigned(3037000493)
    a*old + b
end

function randfast(k, x::Uint)
    x = lcg(x)
    1 + rem(x, k) % Int, x
end

function t4(N, R)
    state = rand(Uint)
    for i = 1:N
        x, state = randfast(R, state)
    end
end

But be careful, if the range is (really) big.

m = div(typemax(Uint),3)*2

julia> mean([rand(1:m)*1.0 for i in 1:10^7])
6.148922790091841e18

julia> m/2
6.148914691236517e18

but (!)

julia> mean([(rand(Uint) % m)*1.0 for i in 1:10^7])
5.123459611164573e18

julia> 5//12*tm
5.124095576030431e18

Upvotes: 2

Nils Gudat
Nils Gudat

Reputation: 13800

Note that in 0.4, int() is deprecated, and you're aske to use round() instead. 

function t2(N,n_mol)
  for i = 1:N
    round(rand()*n_mol)
  end
end

gives 0.27 seconds on my machine (using Julia 0.4).

Upvotes: 1

Related Questions