黄昏より暗きもの、血の流れより赤きもの

読者です 読者をやめる 読者になる 読者になる

黄昏より暗きもの、血の流れより赤きもの

自分の好きな事を好きなように書いて行きます。

【数値計算】ガウス・ルジャンドル法(算術幾何平均と楕円積分)に依る、円周率の近似値の計算【C言語】

C C言語 数値計算

はじめに

今日は円周率の数値計算法のひとつガウスルジャンドル法について説明したい。ガウスルジャンドル法とは算術幾何平均と楕円積分との関係性を利用し、円周率の近似値を得るための方法である。PCのベンチマークソフト「Super Pi(金田康正のホ−ムペ−ジにようこそ! )」にて採用されたアルゴリズムとしても有名である。

前半:円周率計算を使ったPCベンチマークの例

円周率の近似値計算は、パソコンやGPGPU*1や分散処理のベンチマークとして登場する。そこで記事の前半ではベンチマークの例を多数紹介する。ベンチマークが行われた例と、その時使われた式について紹介していく。

例1:Computing Digits of π with CUDA

例として「Computing Digits of π with CUDA」という記事を紹介しよう。使われている計算式「Bellard formula」もこりゃまた複雑。式を見てどういう部分が並列処理に適しているのか皆目検討つかずだ。


f:id:program_study:20160206215048p:plain

ちなみに「four NVIDIA GTX 690 graphics cards」「two NVIDIA GTX 680 graphics cards」「24 computers (at Santa Clara University Design Center) with one NVIDIA GTX 570 graphics card each.」と大学のスパコンと高級なグラフィックボードを刺しまくったときとの結果が比較しているようだ。

例2:Computing PI to 130 000 000 000 decimal digits with y-cruncher..

次に紹介するのはπの計算に「Chudnovsky Formula(wikipedia)」を用いた動画を紹介しよう。まず式を紹介すると、


f:id:program_study:20160206215209p:plain

動画の一部始終は以下のとおり。



Computing PI to 130 000 000 000 decimal digits with y-cruncher..

後半:ガウスルジャンドル法の説明

前半では式1つで円周率を計算する方法を紹介したが、後半ではいよいよガウスルジャンドル法の説明をしようと思う。説明に当たりPDFファイル「算術幾何平均と円周率 上越教育大学 中川 仁」を参考にした。

算術幾何平均AGM(a,b)について

まず2つの正の実数a,b(a>b)がある。このとき相加平均をA(a,b)、相乗平均をG(a,b)するとA(a,b)G(a,b)は以下の式で与えられる。

A(a,b) = \frac{a+b}{2}G(a,b) = \sqrt{ab}

さらにn=0,1,2,…と数列a_nb_nを以下のように定める。

a_{n+1} = \frac{a_n+b_n}{2} b_{n+1} = \sqrt{a_nb_n}

このとき極限値 \lim_{n \to \infty} a_n = αと、 \lim_{n \to \infty} b_n = βが存在する事が知られている。ここで α = \lim_{n \to \infty} a_{n+1} = \lim_{n \to \infty} \frac{a_n+b_n}{2} = \frac{α+β}{2}よりα=βを得る。この極限値αab算術幾何平均といい、AGM(a,b)で表す。

ガウスの公式について

この算術幾何平均と楕円積分というものを経由し、円周率の近似値は以下のように求まる(ガウスの公式)。

各数列をa_0 =1b_0 = c_0 = \frac{1}{\sqrt{2}} a_{n+1} = \frac{a_n+b_n}{2}b_{n+1} = \sqrt{a_nb_n} c_{n+1} = \frac{a_n-b_n}{2} n=0,1,2,…と定める時、円周率の近似値は以下で与えられる。

\pi = \frac{2AGB(1,\frac{1}{\sqrt{2}})}{1 - Σ2^{n-1}{c_n}^2}

ここで分母の値をs_n = {1 - Σ2^{n-1}{c_n}^2}とおく。
a_{n+1} +c_{n+1} = a_n<=>c_{n+1} = a_n - a_{n+1}となるので、
s_{n+1} = s_n - p_n * (a_n - a_{n+1}),p_{n+1} = 2p_nとおくと分母を計算することができる。

次に分子をみると、AGM(a,b)=aまたはbなので、相加相乗平均の式を見返すことにより、(a+b)^2 \geq 4abとなり、最終的に以下を計算すればよい。

\pi = \frac{(a_n+b_n)^2}{4s_n} (\pi = \frac{a_nb_n}{s_n} から変形している。)

まとめ

ここまでの計算法をまとめると以下のようになる。

[1]:数列a_n,b_n,s_n,p_nの初期値を以下に定める。
a_0 =1  b_0= \frac{1}{\sqrt{2}}  t_0= \frac{1}{4}  p_0= 1
[2]:数列a_n,b_n,s_n,p_nを以下に定め、ループによりこの計算をn回繰り返す。
a_{n+1} = \frac{a_n+b_n}{2} b_{n+1} = \sqrt{a_nb_n} s_{n+1} = s_n - p_n * (a_n - a_{n+1}) p_{n+1} = 2p_n
[3]:\pi = \frac{(a_n+b_n)^2}{4s_n}の値を出力する。

プログラム

C言語


最後に

C言語版で反復回数nの回数を3回に指定して出力したところ、πの近似値は3.1415927となった。

参考文献

今回の記事を書く上で、まず「ガウス・ルジャンドルの方法での円周率:きしだのはてな」をみながらプログラムを完成させた。しかしながらその原理の説明が書かれている文献を探すのに時間がかかり、ようやく公開できた事を嬉しく思う。

最後に今回は自分の理解できなかった部分も多数あり、楕円積分の説明及びガウスの公式を求めるまでの手順については説明できなかった。これらの詳しい内容を知りたい方は以下を読んで欲しい。

*1:パソコンにやらせるべき演算をグラフィックボードにやらせよう!という必殺技。Intel Core i7などよりも効果が期待される事がある。