円周率πをパソコンで求める

以前、OpenOfficeのCalcという表計算ソフトを使って、正多角形の辺の長さから円周率を求めてみようという記事を書いた。しかし、結果的に言えば、それはCalcにあらかじめ必要な精度で入っている円周率の値を使ってごまかしているに過ぎない。

今回は、その反省を生かしつつ、もっと数学的に解く方法ふと思いついた（誰でも思いつくだろうが）ので試してみる。

使うのは三角関数の逆関数のテーラー展開。とりあえずArctanを使ってみる。展開できるとかできないとか細かいことはとりあえず置いておく。

(Arctan x)'=1/(1+x^2)=1-x^2+x^4-x^6+x^8-...

※y=Arctan xを考えて、
(Arctan x)'=1/(tan y)'=(cos y)^2=1/(1+(tan y)^2)=1/(1+x^2)

各項を積分して
Arctan x =x-1/3x^3+1/5x^5-...

x=1とすると、
π/4=1-1/3+1/5-1/7+...
π=4(1-1/3+1/5-1/7+...)

ここから先をパソコンにやらせようと思う。

Java Script版もあとで公開しようと思うが、とりあえずCで書いてみた。
各項の±はべき乗でやっても良いんだが、Cだとなんかめんどくさそうなのでk=-1 or 1をおいてある。


#include<stdio.h>

main(){
float pi;
int i,n,k;
i=1;
k=1;
pi=0;
n=101;
while(i<=n){
pi=pi+1.00000000*k/i;
printf("i=%d \t pi=%8.7f \n",i,4*pi);
i=i+2;
k=-1*k;
}
}

このプログラムだと第101項目まで計算しているんだが3.14にはほど遠い。

仕方ないのでこんなずらーっとしたのを見てもうれしくないので、最終結果だけ出るようにいじった上で、記録に挑戦
なお、Wikipedia等によれば正式な値はπ = 3.14159265...だという

i=1001 π=3.1435893
i=100001 π=3.1416154

まだ計算は遅くなったりしていないので大丈夫とみて、さらに項数を増やす。
i=100000001 π=3.1415968

そろそろ結果が出るのにワンテンポ空きが出るようになったので、部屋の暑さも考慮してそろそろやめる。というかfloatでなくてdoubleの方が良いのか・・改良してみよう。

もしかしたらArctanよりも速くπに収束するのがあるかもしれない。それは続編にとっておくことにするが、たった1000ほどで3.14にお目にかかれるとはなかなか意外。

なお、自分で試す場合、nを変えることで精度を変えられる（n→大⇒精度→良）

--精度を上げてみたver.2
double pi,m;

int i,n,k;

i=1;

k=1;

n=1000001;

pi=0.0000000000;

while(i<=n){

m=i;

pi=pi+k/m;

i=i+2;

k=-1*k;

}

printf("i=%d \t pi=%8.7f \n",i-2,4*pi);

有効数字がさらに一桁向上。四捨五入だと思うと表示されている分は完璧。

＊実際に試す方法
Cのコンパイラを入手する。Microsoft Visual StudioのExpress版が良いと思う。(C++)
ファイルを新規作成して.cの拡張子で保存、コンパイルする。
コマンドプロンプトで、プロジェクトのフォルダまで移る。
（設定を変えていなければ、cdを使って、Documents-Visual...-Projects-プロジェクト名）
その後Debugフォルダに移る。

プロジェクト名がsampleとすると、
C:\＊＊＊\sample\Debug>の状態になる。
そこで、sample.exeと入力して[Enter]

n=1000001くらいだとCore2duo 2GHz・省電力モードで1秒ほどかかった。最近のCPUだともっと速いと思うが、調子に乗ってnを増やすと暴走するかもしれないので、自己責任でやってほしい。

（追記）
ちゃんとした解析の本で調べた結果を書いておく。
・今回用いたArctanのテーラー展開でπを求めるのは、収束が遅いため実用的でない。（←確かに・・）
・π/4=4Arctan1/5-Arctan1/239だと比較的速く収束する。（←なぜ？まぁArctanを使うのは誤ってないようだ。）
・現在はこれより速い計算方法がある（らしい←どうやるのか・・・？？）

参考　杉浦光夫「解析入門Ⅰ」東京大学出版会