円周率を求める

概要

2005/06/24

昔のサイトにて、戯言に書かれた円周率算出公式π=lim[θ⇒0](180/θ)sinθは、三角形の面積公式を使った円周率計算にあるような結論となったために、他の展開公式を用いることにした。

独自の関数で実装した円周率計算はこちら

有名なπ展開公式

2005/06/24

ニュートン(速)

π=6(1/2+1/2・3・23+1・3/2・4・5・25+1・3・5/2・4・6・7・27+1・3・5・7/2・4・6・8・9・29…)

グレゴリー・ライプニッツ(遅)

π/4=1-1/3+1/5-1/7+1/9-1/11+1/13+1/15…

マチン(速)

π/4=4arctan(1/5)-arctan(5/239)を展開し

π/4=4{1/5-1/3・53+1/5・55-1/7・57…}-{1/239-1/3・2393+1/5・2395-1/7・2397…}

ブラウンケル(遅)

4/π=1+1/(2+9/(2+25/(2+49/(2+81/(2+121/(2…

ウォーリス(遅)

π=2・(2・2・4・4・6・6・8・8・10・10…)/(1・3・3・5・5・7・7・9・9・11…)

(とりあえず)面倒ではない収束の遅いとされる展開公式が、いかなる性能かを実験してみよう。 ただし、注意してほしいのは面白半分に数値を上げると、計算時間がかかりすぎてフリーズしてしまうこと。

JavaScriptによる円周率簡易計算プロラム

グレゴリー・ライプニッツ

ループ回数

π

20万回のループをやって小数点以下5桁正解。

ウォーリス

ループ回数

π

172回を超えると数値の限界を超えてしまいエラーにorz。(プログラムの計算方式を変えれば解決可能)

とまぁ2つの収束速度の遅い公式を挙げた。当然桁数も限られるので計算不能なんてことも起こる。次にコンピュータによる計算に利用されたマチンの公式。

マチン

ループ回数

π

12回のループで20桁全て合いやがった。流石マチン。

資料

使用された3つのプログラムの内容。

// グレゴリー・ライプニッツ

function GR(CL) {
var R = eval(CL.R.value);
var PT = 0;
for( i=1 ; R>i ; i++ ){
var P1 = Math.pow(-1,i+1)*(1/(2*i-1));
var PT = PT+P1;
var PI = PT*4;
CL.PI.value = PI ;
}}

// ウォーリス

function UO(CL) {
var R = eval(CL.R.value);
var PT1 = 1;
var PT2 = 1;
for( i=2 ; R>i ; i++ ){
var P1 = Math.floor(i/2);
var Pno = P1*2;
var Pun = (P1*2)-1;
var PT1 = PT1*Pno;
var PT2 = PT2*Pun;
var PI = 2*(PT1/PT2)/Pun;
CL.PI.value = PI ;
}}

// マチン

function MC (CL) {
var R = eval(CL.R.value);
var PT1 = 0;
var PT2 = 0;
for( i=1 ; R>i ; i++ ){
var PT1 = PT1+Math.pow(-1,i+1)*(1/((i*2-1)*(Math.pow(5,i*2-1))));
var PT2 = PT2+Math.pow(-1,i+1)*(1/((i*2-1)*(Math.pow(239,i*2-1))));
var PI = 4*(4*PT1-PT2);
CL.PI.value = PI;
}}

ガウス・ルジャンドル

2005/12/15

近年ではマチンなどの逆正接を用いた公式よりFFTを用いた計算方法が用いられる。 代表的なものが以下に示すガウス・ルジャンドルの公式である。ちなみにこれはかの有名なスーパーπでも使われている。 ガウス・ルジャンドルは数学的公式というよりは算法(アルゴリズム)である。

JavaScriptによる計算方法を示す。

function Gauss( max )
{
   var a = 1;
   var b = 1 / Math.sqrt( 2 );
   var t = 1;
   var x = 4;
   for( var i=0 ; i<max ; i++ )
   {
        y = a;
        a = (a + b) / 2;
        b = Math.sqrt(y * b);
        t -= x * Math.pow(y - a, 2);
        x *= 2;
   }
   return Math.pow(a + b, 2) / t;
}

ガウス・ルジャンドル

ループ回数

π

この方法で円周率を延々と求める企画が巨大数演算/HugeIntで行われている。2005/12/15現在JavaScript上で円周率1万桁の計算に成功している。

関連記事

2008/08/04

• 乱数を用いた円周率計算
• 巨大数演算による円周率計算
• 円周率展開公式
• 巨大数演算