多倍長精度計算パッケージ Rmpfr

http://minato.sip21c.org/seek_pi_by_Archimedes.R
中澤さんの，「100ビット実数を使って円周率を10桁正しく出すには正何角形まで必要かを求めるためのRコード」について

Rmpfr というパッケージの使い方がこれだけでわかり，ありがたかった。

もともと，計算時間がかかるというプログラムではなかったが，ちょっとした工夫で計算時間が短くなることがわかったので，記録しておこう。

末尾に追記あり。

中澤さんのオリジナルプログラム

func1 = function() {
    library(Rmpfr)
    n = 6
    p = sqrt(mpfr(3, 100)) / mpfr(2, 100)
    a = sqrt(mpfr(1, 100) - p ^ 2)
    ad = a / p
    L = a * n
    M = ad * n
    while ((M - L) > mpfr(1e-10, 100)) {
        n = n * mpfr(2, 100)
        p = sqrt((mpfr(1, 100) + p) / mpfr(2, 100))
        a = sqrt(mpfr(1, 100) - p ^ mpfr(2, 100))
        ad = a / p
        L = a * n
        M = ad * n
    }
    invisible(c(n, L, M))
}

毎回 mpfr を使うのを避ける

func2 = function() {
    library(Rmpfr)
    n = mpfr(6, 100)
    one = mpfr(1, 100)
    two = mpfr(2, 100)
    epsilon = mpfr(1e-10, 100)
    p = sqrt(mpfr(3, 100)) / two
    a = sqrt(one - p ^ 2)
    ad = a / p
    L = a * n
    M = ad * n
    while ((M - L) > epsilon) {
        n = n * two
        p = sqrt((one + p) / two)
        a = sqrt(one - p ^ two)
        ad = a / p
        L = a * n
        M = ad * n
    }
    invisible(c(n, L, M))
}

混合演算の場合は，自動的に mpfr してくれるようなので，そのまま定数を書く

func3 = function() {
    library(Rmpfr)
    n = 6
    epsilon = mpfr(1e-10, 100)
    p = sqrt(mpfr(3, 100)) / 2
    a = sqrt(1 - p ^ 2)
    ad = a / p
    L = a * n
    M = ad * n
    while ((M - L) > epsilon) {
        n = n * 2
        p = sqrt((1 + p) / 2)
        a = sqrt(1 - p ^ 2)
        ad = a / p
        L = a * n
        M = ad * n
    }
    invisible(c(n, L, M))
}

更に，重箱の隅をつつき，計算時間を搾り取る

func4 = function() {
    library(Rmpfr)
    n = 6
    epsilon = mpfr(1e-10, 100)
    p = sqrt(mpfr(3, 100)) / 2
    a = sqrt(1 - p ^ 2)
    ad = a/p
    L = a * n
    M = ad * n
    while ((M - L) > epsilon) {
        n = n + n
        p = (1 + p) * 0.5
        a = sqrt(1 - p)
        p = sqrt(p)
        L = a * n
        M = L / p
    }
    invisible(c(n, L, M))
}

library(rbenchmark)
benchmark(func1(), func2(), func3(), func4())

結局の所，mpfr が不要な定数はそのまま書くのが速いらしい。

> benchmark(func1(), func2(), func3(), func4())
     test replications elapsed relative user.self sys.self
1 func1()          100  22.560    2.429    22.550    0.146
2 func2()          100  13.763    1.482    13.748    0.075
3 func3()          100  12.017    1.294    12.003    0.066
4 func4()          100   9.287    1.000     9.277    0.051

結果が正しいことを確認

> print(func1())
3 'mpfr' numbers of precision  100   bits
[1]                            786432 3.1415926535814376697324981306951 3.1415926536065043758325431868181
> print(func2())
3 'mpfr' numbers of precision  100   bits
[1]                            786432 3.1415926535814376697324981306951 3.1415926536065043758325431868181
> print(func3())
[[1]]
[1] 786432

[[2]]
'mpfr1' 3.1415926535814376697324981306951

[[3]]
'mpfr1' 3.1415926536065043758325431868181

> print(func4())
[[1]]
[1] 786432

[[2]]
'mpfr1' 3.1415926535814376697324981306951

[[3]]
'mpfr1' 3.1415926536065043758325431868181

Julia で書くと以下のようになる。

function　func_julia(bits=256, epsilon=eps())
    setprecision(BigFloat, bits)
    n = 6
    p = sqrt(big(3)) / 2
    a = sqrt(1 - p ^ 2)
    ad = a/p
    L = a * n
    M = ad * n
    while (M - L) > epsilon
        n = 2n
        p = (1 + p) * 0.5
        a = sqrt(1 - p)
        p = sqrt(p)
        L = a * n
        M = L / p
    end
    (n, L, M)
end

@time n, L, M = func_julia(100, 1e-10)

n # 786432
L # 3.1415926535814376697324981306951
M # 3.1415926536065043758325431868181

func4 に比べると ほぼ 1000 倍速い。