円周率を Machin’s formula で求める

Java

一年近く前に円周率をモンテカルロ法を利用して求めるプログラムを組んだ。

今回は マチンの公式 ( Machin’s formula ) を利用してみようと思う。

そして桁数と精度も上げてみたいと思ってこの方法を選択してみた。

まず、マチンの公式とは次のようになっている。

ifrm0001

この公式の証明は三角関数の公式を利用して証明できます。

倍角公式を二度使い、加法定理でとどめを刺して終わりです。

詳しくはググってみてください。

それではプログラムをこれにそって組んでみましょう。

プログラムを組むにあたって都合の良いように式を変えてあります。( 6 行目 )

このプログラムでは arcTan(double x) の n 番目 (i + i + 1) の項が 0 になるまで計算をします。

十分小さな値をとるまでに変更する場合は、if 文の条件は絶対値をとるように変更すればいいでしょう。

さて、このプログラムの実行結果を確認してみましょう。

3.141592653589794

最後の桁が 3 であるところが 4 になってます。

モンテカルロ法より速くて精度は高いですね。

でも、もっと桁数をとりたくないでしょうか?

そこで、このプログラムを次のように変更してみました。

百万桁の計算をさせるために BigDecimal を使ってみました。

基本的に先ほどのプログラムと同じです。

ちょっと処理時間が気になるので時間計測をいれてます。

と言っても百万桁も計算させると凄く時間がかかってしまうのは解りきっていることなので一万桁で試してみます。

普段使う機会が少ない BigDecimal を使ってますのでちょっとおかしいところがあるかもしれませんが気にしないでね。(^_^;)

あと、桁数を多くとるので 0 収束を使うと非常に時間がかかってしまうので誤差がでないように適当にループさせてます。

それではいつものように未熟で汚いコードですがご覧ください。

少しでも計算時間が短くなればいいなぁって思って悪あがきしてます。( ̄。 ̄;)

基本的な考え方は double 仕様のものと一緒で悪あがきの並行処理を取り入れただけです。

実行結果の抜粋を載せておきます。

<— 円周率計算開始 —>

10000桁の円周率は、
3.1415926535897932384626433832795028841971693993751058209749445923078164062

4988727584610126483699989225695968815920560010165525637568

10000桁の円周率計算時間は、0時間3分22秒4562651430000244

<— 円周率計算終了 —>

ちゃんと計算できてますね。

一番最後は丸めによって四捨五入されてます。

もっと効率よくプログラムを組むには収束値を設定したほうがいいのかなぁ・・・?

でも、速くても計算が甘くなるのはいやだ!

とりあえず昨年作ったモンテカルロ法で円周率を求めるプログラムより遙かに良いのは確かだ!

あと気のせいかもしれないが BigDecimal の除算って遅く感じる。

ゴールデンウィークの最後の休日にこんなことしているってのは平和でいいね。(^_^)

平成25年5月8日(水) 追記 計算のループの回数を決めている LIMIT_NUMBER の値を変更しました。

10000 桁で確認したところ問題なさそうでした。

ただし、それより大きな桁は確認していませんのであしからず。(;´Д`)

Hatena タグ:
このエントリーをはてなブックマークに追加
Share on Facebook

« »

Comment

  1. Good post. I’m going through some of these issues as well..

    返信

Leave a Reply

* が付いている項目は、必須項目です!

*

次の HTML タグと属性を利用できます: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code class="" title="" data-url=""> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong> <pre class="" title="" data-url=""> <span class="" title="" data-url="">

Trackback URL