在 Julia 中生成范围内随机整数的有效方法

Question

我正在做 MC 模拟，我需要在 1 和变量上限 n_mol 之间的范围内生成随机整数

执行此操作的特定 Julia 函数是 rand(1:n_mol)，其中 n_mol 是一个整数，随着每次 MC 迭代而变化。问题是这样做很慢......（可能是 Julia 开发人员需要解决的问题）。所以，我没有使用那个特定的函数调用，而是考虑在 [0,1) 中生成一个随机浮点数乘以n_mol，然后得到结果的整数部分：int(rand()*n_mol) 现在的问题是int()四舍五入，所以我可以得到0和n_mol之间的数字，我不能得到0...所以我目前使用的解决方案是使用ifloor并添加@987654333 @，ifloor(rand()*n_mol)+1，比第一个快很多，但比第二个慢。

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)

那么，有什么方法可以在第二次测试附近的速度下更快地做到这一点？ 这很重要，因为 MC 模拟进行了超过 1e10 次迭代。 结果必须是整数，因为它将用作数组的索引。

Answer 1

鉴于以下两个考虑，rand(r::Range) 代码非常快。首先，julia 调用 52 位 rng 两次以获得随机整数，调用 52 位 rng 一次以获得随机浮点数，这给出了一些记账因素 2.5。第二件事是

(rand(Uint) % k) 

仅在 0 到 k-1 之间均匀分布，如果 k 是 2 的幂。这是通过拒绝抽样来处理的，这或多或少地解释了剩余的额外成本。

如果速度非常重要，您可以像 Julia 一样使用更简单的随机数生成器并忽略这些问题。例如没有拒绝采样的线性同余生成器

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

但要小心，如果范围（真的）很大。

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

但是（！）

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

julia> 5//12*tm
5.124095576030431e18

Answer 2

请注意，在 0.4 中，int() 已弃用，您需要改用 round()。

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

在我的机器上提供 0.27 秒（使用 Julia 0.4）。