Julia でビュフォンの針

以前に掲載したが，π を使わないシミュレーション

n = 100000000 で，

Python 35.806 秒
R      14.543 秒
Julia   5.836 秒

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

★★　Python

import numpy as np
def buffon(n = 100000000): # 取りあえずの試行回数（実際はこれより少なくなる）
    d = 10 # 平行線の間隔
    l = d / 2 # 針の長さ
    # 以下 5 行で，0 <= theta < π/2 の一様乱数をもとめ，その sine.theta を生成する
    xy = np.random.rand(2, n) # 一辺 1 の正方形内の一様乱数による座標点
    r = (xy**2).sum(axis=0) # 原点からの距離
    xy = xy[:, r <= 1] # 半径 1 の円の内部にある点に限る
    r = r[r <= 1] # 半径 1 の円の内部にある点に限る
    sine_theta = xy[1,] / np.sqrt(r) # 三角関数の定義により
    n = len(r) # sine.theta の個数
    y = np.random.rand(n)*l # 同じ個数の[0, l) の一様乱数（針の中心から，平行線までの距離）
    return 2 * l / d / np.mean(y <= l / 2 * sine_theta) # 条件を満たすものの逆数が π の推定値

from time import time
s = time(); print(buffon(100000000), time() - s) # 3.1409810174102484 35.805713176727295

★★　R

buffon = function(n = 1e+06) { # 取りあえずの試行回数（実際はこれより少なくなる）
 d = 10 # 平行線の間隔
 l = d/2 # 針の長さ
 # 以下 5 行で，0 <= theta < π/2 の一様乱数をもとめ，その sine.theta を生成する
 xy = matrix(runif(2 * n, 0, 1), 2) # 一辺 1 の正方形内の一様乱数による座標点
 r = colSums(xy^2) # 原点からの距離
 xy = xy[, r <= 1] # 半径 1 の円の内部にある点に限る
 r = r[r <= 1] # 半径 1 の円の内部にある点に限る
 sine.theta = xy[2, ]/sqrt(r) # 三角関数の定義により
 n = length(r) # sine.theta の個数
 y = runif(n, 0, l) # 同じ個数の[0, l) の一様乱数（針の中心から，平行線までの距離）
 2 * l/d/mean(y <= l/2 * sine.theta) # 条件を満たすものの逆数が π の推定値
}

system.time(print(buffon(100000000)))
# [1] 3.141178
#    ユーザ   システム       経過  
#     14.543      2.852     17.349 

★★　Julia

function buffon(n = 1000000) # 取りあえずの試行回数（実際はこれより少なくなる）
    # 以下 5 行で，0 <= theta < π/2 の一様乱数をもとめ，その sinθ を生成する
    xy = rand(2, n);  # 一辺 1 の正方形内の一様乱数による座標点
    r = vec(sum(xy .^ 2, dims=1)); # 原点からの距離
    xy = xy[:, r .<= 1]; # 半径 1 の円の内部にある点に限る
    r = r[r .<= 1]; # 半径 1 の円の内部にある点に限る
    sinθ = xy[2, :] ./ sqrt.(r); # 三角関数の定義により
    n = length(r) # sinθ の個数
    n / sum(rand(n) .<= (0.5 .* sinθ)) # 条件を満たすものの逆数が π の推定値
end

@time buffon(100000000) #   5.836088 seconds (30 allocations: 7.269 GiB, 6.80% gc time)
# 3.141511691253939