円周率を計算する その１

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

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

ガウス＝ルジャンドルのアルゴリズム
https://ja.wikipedia.org/wiki/%E3%82%AC%E3%82%A6%E3%82%B9%EF%BC%9D%E3%83%AB%E3%82%B8%E3%83%A3%E3%83%B3%E3%83%89%E3%83%AB%E3%81%AE%E3%82%A2%E3%83%AB%E3%82%B4%E3%83%AA%E3%82%BA%E3%83%A0

なるほど、
初期値

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

反復式

 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を計算

  \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;
}

結果

Start..
 pi = 3.140579250522169
End..

実際は　3.141592... なので
3.14 まではあっている。

以外に簡単な計算式(アルゴリズム)なのだな。
（内容の理解はともかく）

回数を増やせば精度が上がるので、次回は繰り返して計算できるようにソースを改変してみる。