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

円周率を計算する その1

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

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

ガウス=ルジャンドルのアルゴリズム - Wikipedia

なるほど、 初期値

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

    //1回目
    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 まではあっている。

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

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