TIM Labs

円周率を1万桁まで求めてみた

| コメント(0) | トラックバック(0)
さて、円周率を連分数展開して計算することで有理数としては、いくらでも精度をあげられるらしいことが判明した。

しかし、有理数(分数)のままでは困る。10進数に直したい。

ということで、小数部分を、10**n を掛けることで、小数点の位置をずらしてしまおう。
そうして、整数部分だけをtruncateで取り出せば良いはずだ。

ということで、とりあえず小数点以下50桁まで求めてみよう。
そのために、 10 ** 50 を掛けてみる。

Prelude Data.List Data.Ratio> 10**50
1.0e50

想定外の指数表示になってしまった。
有理数にして、きちんと確認してみよう。

Prelude Data.List Data.Ratio> toRational (10**50)
100000000000000007629769841091887003294964970946560 % 1

と思ったら、なんと、1の後ろに0が50個並ばなかった。
雑な近似が行われていることが判明した。
** を使った場合、浮動小数点数としての計算で誤魔化されるみたいだ。

ベキ乗の演算子は、** と ^ の2つがあった。
それでは、^ を使うと、どうなるだろうか。

Prelude Data.List Data.Ratio> 10^50
100000000000000000000000000000000000000000000000000

こちらは、誤差がでないので、こちらを使おう。

どうやら、これで精度は問題ないみたいだ。
連分数を100項まで計算し、10^50との積の整数部分を取り出してみた。

Prelude Data.List Data.Ratio> truncate ((contfracPi' 100) * 10^50)
314159265358979323846264338327950288419716939937510

円周率の値はネット上にゴロゴロしているので比較してみたが、大丈夫のようだ。
ということで、調子に乗って、1000桁まで求めてみた。連分数は1500項までの近似を使った。

Prelude Data.List Data.Ratio> truncate ((contfracPi' 1500) * 10^1000)
3
14159265358979323846264338327950288419716939937510
58209749445923078164062862089986280348253421170679
82148086513282306647093844609550582231725359408128
48111745028410270193852110555964462294895493038196
44288109756659334461284756482337867831652712019091
45648566923460348610454326648213393607260249141273
72458700660631558817488152092096282925409171536436
78925903600113305305488204665213841469519415116094
33057270365759591953092186117381932611793105118548
07446237996274956735188575272489122793818301194912
98336733624406566430860213949463952247371907021798
60943702770539217176293176752384674818467669405132
00056812714526356082778577134275778960917363717872
14684409012249534301465495853710507922796892589235
42019956112129021960864034418159813629774771309960
51870721134999999837297804995105973173281609631859
50244594553469083026425223082533446850352619311881
71010003137838752886587533208381420617177669147303
59825349042875546873115956286388235378759375195778
18577805321712268066130019278766111959092164201989

1000桁を1秒程度で計算してくれるようだ。なかなか収束が良い連分数のようだ。

さらに調子に乗って、10000桁を求めてみた。

Prelude Data.List Data.Ratio> truncate ((contfracPi' 15000) * 10^10000)

しかし、なかなか反応がないので、食事に行ったら、結果が出ていた。

3.
14159265358979323846264338327950288419716939937510
14159265358979323846264338327950288419716939937510
58209749445923078164062862089986280348253421170679
82148086513282306647093844609550582231725359408128

................ 中略 ................

26456001623742880210927645793106579229552498872758
46101264836999892256959688159205600101655256375678

実際は、10001桁の数字が延々と続いていただけだが。

結局、40分くらいかかったようだ。

円周率を追いかけ回すのは、このあたりでお終いにしよう。

トラックバック(0)

トラックバックURL: http://labs.timedia.co.jp/mt/mt-tb.cgi/339

コメントする

このブログ記事について

このページは、fujiが2012年12月 6日 10:00に書いたブログ記事です。

ひとつ前のブログ記事は「フィボナッチ数列の一般項をもう少しだけ速く計算する(ライブラリ活用編)」です。

次のブログ記事は「Haskellでお絵描き」です。

最近のコンテンツはインデックスページで見られます。過去に書かれたものはアーカイブのページで見られます。