円周率を1万桁まで求めてみた
2012年 12月 06日
さて、円周率を連分数展開して計算することで有理数としては、いくらでも精度をあげられるらしいことが判明した。
しかし、有理数（分数）のままでは困る。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個並ばなかった。
雑な近似が行われていることが判明した。
** を使った場合、浮動小数点数としての計算で誤魔化されるみたいだ。
ベキ乗の演算子は、** と ^ の２つがあった。
それでは、^ を使うと、どうなるだろうか。
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
185778053217122680661300192787661119590921642019891000桁を１秒程度で計算してくれるようだ。なかなか収束が良い連分数のようだ。
さらに調子に乗って、10000桁を求めてみた。
Prelude Data.List Data.Ratio> truncate ((contfracPi' 15000) * 10^10000)しかし、なかなか反応がないので、食事に行ったら、結果が出ていた。
3.14159265358979323846264338327950288419716939937510
14159265358979323846264338327950288419716939937510
58209749445923078164062862089986280348253421170679
82148086513282306647093844609550582231725359408128
................ 中略 ................
26456001623742880210927645793106579229552498872758
46101264836999892256959688159205600101655256375678実際は、10001桁の数字が延々と続いていただけだが。
結局、40分くらいかかったようだ。
円周率を追いかけ回すのは、このあたりでお終いにしよう。
＃人気の記事
＃タグ
