円周率を計算する　その１ - なんとなくのブログ

円周率をPCで計算する。実際、円周率計算ソフトなどをつかって円周率を求めることができるが
スクラッチ状態から円周率プログラムを組んだことがなかったので実際にやってみた。

どうやって円周率を求めるのか、そのアルゴリズムを検索。
「円周率　アルゴリズム」でググってみると最初に出てくるのが

なるほど、初期値

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

反復式

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

を何回か繰り返して

PIを計算

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

する。

ループなしで一回のみ計算(n+1回目のみ)してみて円周率をもとめてみた。
以下がそのソース。

#include <math.h>
#include <iostream>
#include <iomanip>

void gauss_Legendre_algorithm()
{
    std::cout << "Start.." << std::endl;

    //初期値
    double a0 = 1.0;
    double b0 = 1.0 / sqrt(2.0);
    double t0 = 1.0 / 4.0;
    double p0 = 1.0;

    //１回目
    double a1 = (a0 + b0) / 2.0;
    double b1 = sqrt(a0 * b0);
    double t1 = t0 - p0 * pow((a0 - a1), 2.0);
    double p1 = 2 * p0;

    //PI計算
    double pi = 0.0;
    pi = pow((a1 + b1), 2.0) / (4.0 * t1);
    std::cout << " pi = "
              << std::fixed
              << std::setprecision(15)
              << pi
              << std::endl;

    std::cout << "End.." << std::endl;
    std::cin.ignore();
}

int main(int argc, char *argv[])
{
    gauss_Legendre_algorithm();
    return 0;
}