Эффективного способ генерации случайных чисел в пределах диапазона в Джулии

Question

Я делаю MC моделирование и мне нужно для генерации случайных чисел в пределах диапазона между 1 и переменным верхним пределом n_mol

Специфической функцией Julia для этого является rand(1:n_mol) где n_mol - целое число, которое изменяется с каждой итерацией MC. Проблема в том, что делать это медленно ... (возможно, проблема для разработчиков Julia). Поэтому вместо того, чтобы использовать этот вызов функции, я подумал о создании случайного поплавка в [0,1], умножьте его на n_mol, а затем получите целочисленную часть результата: int(rand()*n_mol) проблема в том, что int() округляет, чтобы я мог закончить с 0 и , и я не могу получить 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 один раз, чтобы получить случайные поплавки, что дает с некоторой книгой, сохраняющей фактор 2.5. Второе то, что

(rand(Uint) % k) 

только равномерно распределено между 0 до K-1, если к является степенью 2. Это заботится с отбором проб отторжения, это объясняет более или менее оставшиеся дополнительные расходы.

Если скорость чрезвычайно важна, вы можете использовать более простой генератор случайных чисел в качестве Джулии и игнорировать эти проблемы. Например, с линейной конгруэнтной генератор без отказа выборки

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).