「AGM(算術幾何平均)を用いる方法」 で紹介した公式で計算してみよう。 この公式はどうしてこれで円周率が計算できるかを理解するのが難しいが、 アルゴリズムそのものは比較的単純で、 短いプログラムで計算が実行可能である。
ここでは、 鎌田伊織・吉本清夏「 -- 計算法の変遷 --」, 2005年度卒業研究 から拝借したプログラムを紹介する。
GaussLegendre.BAS |
REM GaussLegendre.BAS --- Gauss-Legendre公式によるπの計算 REM 鎌田伊織・吉本清夏,「π-- 計算法の変遷 ---」, 2005年度卒業研究 OPTION ARITHMETIC DECIMAL_high OPTION BASE 0 DIM A(100),B(100),T(100) LET A(0)=1 LET B(0)=1/SQR(2) LET T(0)=1/4 LET MAXN=10 FOR n=1 TO maxn LET a(n)=(a(n-1)+b(n-1))/2 LET b(n)=SQR(a(n-1)*b(n-1)) LET t(n)=t(n-1)-2^(n-1)*(a(n)-a(n-1))^2 PRINT USING "### -%.####^^^^^^":n,PI-(A(n)+B(n))^2/(4*t(n)) NEXT n LET n=maxn LET k=(a(n)+b(n))^2/(4*t(n)) PRINT k END |
GaussLegendre.TXT |
1 1.0134E-0003 2 7.3763E-0009 3 1.8313E-0019 4 5.4721E-0041 5 2.4061E-0084 6 2.3086E-0171 7 1.0586E-0345 8 1.1110E-0694 9 -1.5022E-1000 10 -1.5022E-1000 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420199 |
わずか 回未満の反復で 桁の精度を達成していることに注目。 合っている桁数が反復1回ごとにほぼ2倍になるという、 いわゆる 次の収束をしている (これは例えば、 非線形方程式の解法として紹介した Newton 法の性質と同じ) ので、 あっと言う間に 桁に達成!、と言えるでしょうか。
( の級数展開と比べると、もの凄い違いですね。)