【初学者】Juliaを使って円周率を求める

新しめの言語をお勉強しないと会社でゴミ扱いされそうなので、Juliaを使って簡単な何かを作ってみます。

Table of Contents

Juliaとは

Juliaとは、セクシー女優です。

違います。科学技術計算分野でのプログラミング言語の座をPython, R, MATLABから奪い取ろうとしてる言語です。

色々な会話でちやほやされている言語ですが、お前何がすごいねんを少し掘り下げます。

私動的プログラミング言語なのよ

Juliaは動的プログラミング言語です。コンパイルは実行時に処理されるということです。

なんや、Pythonと同じやんけと思ってましたがところがどっこい。

julialang.orgに次のように書かれていました。

Julia is fast!

Julia was designed from the beginning for high performance. Julia programs compile to efficient native code for multiple platforms via LLVM.

なんでそんなに早いのかというと実行時コンパイル(JIT)だからです。

科学技術計算に多いループ処理を効率よく実行するわけですね。

いやいや、じゃあPythonのPyPy使えばええんやい？って声が聞こえてくるかもしれませんが、ポイントはJITの方式です。

PyPyはJavaでいうところのHotspotみたいな方式(TracingJITというらしいです)でJITを行っているらしく、起動時プロファイルを取りながらHotSpotに対してコンパイルや最適化を実施していくようです。

その点JuliaはLLVMという技術を使って型推論と最適化を行ったのち実行されます。

型が同定できることで強力な最適化がなされるので早いらしいです。

どうでもいいですがLLVMのロゴのドラゴンかっこいいですね。遊戯王カードに出てきそうですね。

遊戯王カード

PythonのNambaもLLVMでJITしてます。なのでもちろん型の指定もできるのでこの手の最適化は強いはず。（早いとは聞いたことがありますが試したことはないです。）

言語としてネイティブにLLVM対応しているのがJuliaのすごいところだと思います。

円周率を計算してみる

別にベンチマークをとってやろうとかそう言うわけではないのですが、取り急ぎ練習として円周率を計算してみたいと思います。

この手の計算には有名どころがいくつかあります。

ライプニッツの公式

超有名どころですね。

summation symbolで書くと

これをJuliaで書くと

function leibniz(n::BigInt, debug=false)
    ret::BigFloat = 0.0
    for i in 0:n
        ret += (-1)^i / (2i + 1)
    end
    return 4ret
end

BigInt, BigFloatが使えるのでこの手の演算は気持ちいいですねー。

julia> leibniz(BigInt(10000)

3.141692643590543213460768320877940222544825752138710733999805489190209879980251

実際に演算してみると、10000ループでも全く精度出ていないことがわかります。

3.14169......

ライプニッツの公式は実際の円周率の計算に用いるには収束が非常に遅いために全く適していないので別の方法も試してみます。

ガウス＝ルジャンドルのアルゴリズム

というのがあるらしいです。

スーパーコンピュータの検収時このアルゴリズムを使って、πの超高精度計算が用いられることが有名です。

初期値

を定義したうえでイテレーション毎に

という計算をすると

とπに近似するという感じです。

これをJuliaで書くと

function gauss_legendre(n::BigInt, debug=false)
    local a::BigFloat = 1.0
    local b::BigFloat = 1.0 / sqrt(2.0)
    local t::BigFloat = 1.0 / 4
    local p::BigFloat = 1.0
    local tmp::BigFloat = 0

    for i in 0:n
        tmp = a
        a = (tmp + b) / 2
        b = sqrt(tmp * b)
        t = t - (p * (a - tmp) * (a - tmp))
        p = 2p
    end
    return (a + b) * (a + b) / 4t
end

おおー直感的ですね。

julia> gauss_legendre(BigInt(10000)

3.14159265358979311120035621508437779160333415473350005136581079848894037301245

私が覚えている範囲3.1415926535...

まではあっていますね。

結論

Juliaをもっと書こうかと思いました。あと、ベンチマークは測ってみます。