理系の憧れ！円周率の求め方　～ガウス＝ルジャンドルのアルゴリズム編～

Tweet

はじめに

前回はライプニッツの公式編を書きました。しかし、それでは収束が遅く、精度がなかなか出にくいということがわかりました。そこで、さらに早いアルゴリズムを用い、どんどん計算していこうと思います。今回は、ガウス様、ルジャンドル様のアルゴリズムをプログラム化していきたいと思います。

環境

OS:Windows 7

IDE:Visual Studio 2015 C++

CPU:Intel Core i7-2600K 3.40GHz

RAM:8.00GB

ちなみに、これはサブPCです。。。

ガウス＝ルジャンドルのアルゴリズムについて

これは漸化式を使ってπを近似する方法です。以下、Wikiより引用。

初期値の設定

$$a_0=1 \quad b_0=\frac {1}{\sqrt{2}} \quad t_0=\frac{1}{4} \quad p_0=1$$

反復式

$$a_{n+1}= \frac{a_n + b_n}{2}\\
b_{n+1} = \sqrt{a_n b_n}\\
t_{n+1} = t_n - p_n(a_n - a_{n+1})^2\\
p_{n+1} = 2p_n$$

πの算出

円周率πは、a,b,tを用いて以下のように近似される。

$$\pi \approx \frac{(a+b)^2}{4t}$$

プログラムにする

まあまあきれいな漸化式なので、普通にプログラム化してしまいます。long doubleでやってみます。最近のパソコンはメモリがたくさん載っているので、floatよりもdoubleを使ったほうがいろいろ効率が良いという話をどこかできいたような…

#include <math.h>
#include <iostream>
#include <iomanip>
using namespace std;

int main(void) {
	long double a, b, t;
	int p, n;
	a = 1.0;
	b = 1.0 / powl(2, 0.5);
	t = 1.0 / 4.0;
	p = 1;
	n = 0;
	cin >> n;
	for (int i = 0; i < n; i++) {
		long double a2, b2, t2;
		int p2;
		a2 = (a + b) / 2.0;
		b2 = powl(a*b, 0.5);
		t2 = t - p*powl(a - a2, 2.0);
		p2 = 2 * p;
		a = a2;
		b = b2;
		t = t2;
		p = p2;
	}
	cout << setprecision(LDBL_DIG) << powl(a + b, 2.0) / (4.0 * t) << endl;
	return 0;
}

単純に書けてしまいました。それでも、ライプニッツの公式と比べると、とてつもなく早くlong doubleの限界まで達してしまいます。wikiによると、オーバークロックマニアの中では有名な「スーパーπ」というソフトにも、このアルゴリズムが使われているそうです。

収束

n=3で、すぐにlong doubleの限界に達してしまいました。

ライブラリの力を借りて、さらに精度を上げていきたいところですね…
近いうちにやりたいと思います。

ありがとうございました。