Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

faster randn by separating out unlikely branch in a function #9132

Merged
merged 1 commit into from
Nov 25, 2014
Merged

faster randn by separating out unlikely branch in a function #9132

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
40 changes: 21 additions & 19 deletions base/random.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,7 @@ type Close1Open2 <: FloatInterval end

@inline rand_ui52_raw_inbounds(r::MersenneTwister) = reinterpret(UInt64, rand_inbounds(r, Close1Open2))
@inline rand_ui52_raw(r::MersenneTwister) = (reserve_1(r); rand_ui52_raw_inbounds(r))
@inline rand_ui52(r::MersenneTwister) = rand_ui52_raw(r) & 0x000fffffffffffff
@inline rand_ui2x52_raw(r::MersenneTwister) = rand_ui52_raw(r) % UInt128 << 64 | rand_ui52_raw(r)

function srand(r::MersenneTwister, seed::Vector{UInt32})
Expand Down Expand Up @@ -943,36 +944,37 @@ const ziggurat_nor_r = 3.6541528853610087963519472518
const ziggurat_nor_inv_r = inv(ziggurat_nor_r)
const ziggurat_exp_r = 7.6971174701310497140446280481

@inline randi(rng::MersenneTwister=GLOBAL_RNG) = reinterpret(Uint64, rand(rng, Close1Open2)) & 0x000fffffffffffff

function randmtzig_randn(rng::MersenneTwister=GLOBAL_RNG)
@inbounds begin
r = rand_ui52(rng)
rabs = int64(r>>1) # One bit for the sign
idx = rabs & 0xFF
x = ifelse(r % Bool, -rabs, rabs)*wi[idx+1]
rabs < ki[idx+1] && return x # 99.3% of the time we return here 1st try
return randmtzig_randn_unlikely(rng, idx, rabs, x)
end
end

# this unlikely branch is put in a separate function for better efficiency
function randmtzig_randn_unlikely(rng, idx, rabs, x)
@inbounds if idx == 0
while true
r = randi(rng)
rabs = int64(r>>1) # One bit for the sign
idx = rabs & 0xFF
x = ifelse(r % Bool, -rabs, rabs)*wi[idx+1]
if rabs < ki[idx+1]
return x # 99.3% of the time we return here 1st try
elseif idx == 0
while true
xx = -ziggurat_nor_inv_r*log(rand(rng))
yy = -log(rand(rng))
if yy+yy > xx*xx
return (rabs & 0x100) != 0x000000000 ? -ziggurat_nor_r-xx : ziggurat_nor_r+xx
end
end
elseif (fi[idx] - fi[idx+1])*rand(rng) + fi[idx+1] < exp(-0.5*x*x)
return x # return from the triangular area
end
xx = -ziggurat_nor_inv_r*log(rand(rng))
yy = -log(rand(rng))
yy+yy > xx*xx && return (rabs >> 8) % Bool ? -ziggurat_nor_r-xx : ziggurat_nor_r+xx
end
elseif (fi[idx] - fi[idx+1])*rand(rng) + fi[idx+1] < exp(-0.5*x*x)
return x # return from the triangular area
else
return randmtzig_randn(rng)
end
end

function randmtzig_exprnd(rng::MersenneTwister=GLOBAL_RNG)
@inbounds begin
while true
ri = randi(rng)
ri = rand_ui52(rng)
idx = ri & 0xFF
x = ri*we[idx+1]
if ri < ke[idx+1]
Expand Down