裏 RjpWiki

文字通り,RjpWiki の裏を行きます
R プログラム コンピュータ・サイエンス 統計学

広告

※このエリアは、60日間投稿が無い場合に表示されます。記事を投稿すると、表示されなくなります。

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

2014年07月30日 | ブログラミング

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

ジャンル:
ウェブログ
コメント   この記事についてブログを書く
« Excel だって,賢く使おうョ | トップ | 3 次元空間の単位ベクトルを... »
最近の画像もっと見る

コメントを投稿

ブログ作成者から承認されるまでコメントは反映されません。

ブログラミング」カテゴリの最新記事

関連するみんなの記事