円周率を計算する その２

前回の続き

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

前回は一回しか計算できないプログラムなので 今回は　n回　計算できるようにソースを改変 以下がそのソース。

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

class Parameter
{
  public:
    double a;
    double b;
    double t;
    double p;

    Parameter() : a(1.0), b(1.0 / sqrt(2.0)), t(1.0 / 4.0), p(1.0)
    {
    }
};

double gauss_Legendre_algorithm(Parameter &param)
{
    double a = (param.a + param.b) / 2.0;
    double b = sqrt(param.a * param.b);
    double t = param.t - param.p * pow((param.a - a), 2.0);
    double p = 2.0 * param.p;

    param.a = a;
    param.b = b;
    param.t = t;
    param.p = p;

    //pi計算
    return pow((a + b), 2.0) / (4.0 * t);
}

int main(int argc, char *argv[])
{
    std::cout << "Start.." << std::endl;
    Parameter param;
    double pi;

    for (int i = 1; i < 6; i++)
    {
        pi = gauss_Legendre_algorithm(param);
        std::cout << "pi = "
                  << std::fixed
                  << std::setprecision(15)
                  << pi
                  << std::endl;
    }

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

結果

Start..
pi = 3.140579250522169
pi = 3.141592646213543
pi = 3.141592653589794
pi = 3.141592653589794
pi = 3.141592653589794
End..

実際は
3.141592653589793238462643... なので
3.14159265358979 まではあっている。

しかし！！　５回計算してみたのだが３回目からはすべて同じ答えになっていて
これ以上続けても無意味！！

doubleの仮数部は52bitなので
double 型の精度（有効桁数）は２進数にして 53 (=52+1) 桁であり，１０進数では約 15 桁となって
これ以上の精度は出せない。

う〜ん、意外に長大な円周率は簡単には表現できないらしい。
ではどうする？

そもそも double 型の精度 がよろしくないのであって、無制限の精度で計算できればいいのだ。
多倍長doubleを自分で組むのは最初から諦めて（世のすごい人たちの成果を利用させてもらう）

多倍長ライブラリ GNU Multi-Precision Library（GMP)
The GNU MP Bignum Library　を利用、次回は精度を上げてみる。