JuliaLangAC2014 -- Day 18 sample codes

collatz.jl

function collatz_branch(N::Integer)
  i = 0
  while N > 1
    N = iseven(N) ? div(N,2) : 3N+1
    i += 1
  end
  return i
end

function collatz_select(N::Integer)
  i = 0
  while N > 1
    N = ifelse(iseven(N), div(N,2), 3N+1)
    i += 1
  end
  return i
end

test_branch(Ns::Vector{Int}) =  map(collatz_branch, Ns)
test_select(Ns::Vector{Int}) =  map(collatz_select, Ns)

gc_disable()

test_branch([1:10])
test_select([1:10])

const N = (length(ARGS) > 0 ? int(ARGS[1]) : 10000)
const Ns = [1:N]

@time test_branch(Ns)
@time test_select(Ns)

random.jl

## Base.randn() がどうやって高速化されたか

## 必要なパラメータや関数を持ってくる
import Base.Random: GLOBAL_RNG, rand_ui52
import Base.Random: ki, wi, fi, ke, we, fe
import Base.Random: ziggurat_nor_r, ziggurat_nor_inv_r, ziggurat_exp_r

## 初期バージョン
@inline function randn_first(rng::AbstractRNG=GLOBAL_RNG)
    @inbounds begin
        r = rand_ui52(rng)
        rabs = int64(r>>1)
        idx = rabs & 0xFF

        # Point 1
        # 分岐命令
        x = (r % Bool ? -rabs : rabs)*wi[idx+1]

        rabs < ki[idx+1] && return x # 99.3% の確率で、ここでreturn

        # Point 2
        # ここに来る確率は 0.7 %
        @inbounds if idx == 0
          while true
            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 randn_first(rng)
        end
    end
end


## Point 1:
## 分岐命令(if ... end / ?: ) をifelse() 関数に置き換え

@inline function randn_ifelse(rng::AbstractRNG=GLOBAL_RNG)
    @inbounds begin
        r = rand_ui52(rng)
        rabs = int64(r>>1)
        idx = rabs & 0xFF

        # Point 1
        # ifelse() 関数（分岐命令なしで値が選択される！）
        x = ifelse(r % Bool, -rabs, rabs)*wi[idx+1]

        rabs < ki[idx+1] && return x # 99.3% の確率で、ここでreturn

        @inbounds if idx == 0
          while true
            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 randn_ifelse(rng)
        end
    end
end

## Point 2:
## およそ0.7% でしか到達しないのにやたら長い部分を
## 別関数に分離

@inline function randn_separate(rng::AbstractRNG=GLOBAL_RNG)
    @inbounds begin
        r = rand_ui52(rng)
        rabs = int64(r>>1)
        idx = rabs & 0xFF
        x = (r % Bool ? -rabs : rabs)*wi[idx+1]
        rabs < ki[idx+1] && return x # 99.3% の確率で、ここでreturn

        # Point 2
        # ここに来る確率は0.7 %
        # 関数呼び出しに置き換えて関数の長さを減らす
        return separate_unlikely(rng, idx, rabs, x)
    end
end

# Point 2
# たまにしか使われないくせに長い部分
function separate_unlikely(rng, idx, rabs, x)
    @inbounds if idx == 0
        while true
            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 randn_separate(rng)
    end
end


## 20141218 現在のrandn()
## Point 1+2
@inline function randn_final(rng::AbstractRNG=GLOBAL_RNG)
    @inbounds begin
        r = rand_ui52(rng)
        rabs = int64(r>>1)
        idx = rabs & 0xFF

        # Point 1
        x = ifelse(r % Bool, -rabs, rabs)*wi[idx+1]

        rabs < ki[idx+1] && return x # 99.3% の確率で、ここでreturn

        # Point 2
        return randn_unlikely(rng, idx, rabs, x)
    end
end

function randn_unlikely(rng, idx, rabs, x)
    @inbounds if idx == 0
        while true
            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 randn_final(rng)
    end
end

gc_disable()


randn_first()
randn_ifelse()
randn_separate()
randn_final()

const N = (length(ARGS) > 0 ? int(ARGS[1]) : 10000 * 10000)

println("初期バージョン")
@time for i in 1:N
  randn_first()
end

println()
println("Point 1: ifelse() 関数を使うことで分岐命令を消す")
@time for i in 1:N
  randn_ifelse()
end

println()
println("Point 2: めったに辿り着かないくせに長い分岐を関数呼び出しに置き換え")
@time for i in 1:N
  randn_separate()
end

println()
println("Point 1+2: 20141218 現在最新バージョン")
@time for i in 1:N
  randn_final()
end