random - 在 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 并添加一个 1 ， 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 一次来获取随机 float ，这给出了一些账簿保留的因子 2.5。第二件事是

(rand(Uint) % k) 

如果 k 是 2 的幂，

仅均匀分布在 0 到 k-1 之间。这是通过拒绝采样来解决的，这或多或少解释了剩余的额外成本。

如果速度极其重要，您可以使用 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