円周率を計算する その３

前回の続き

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

前回はdouble型の精度の問題で15桁までしか計算できなかった。 今回は多倍長ライブラリ(GMP)を使って任意の精度まで計算するように変更。
以下がそのソース。
digits　変数で任意に精度を設定　今回は 200 でやってみる。

//★精度 ここを調整する
int digits = 200;

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

mpf_class power(mpf_class x, int n)
{
    mpf_class powval(1);
    for (int num = n; num != 0; --num)
    {
        powval *= x;
    }
    return powval;
}

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

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

mpf_class gauss_Legendre_algorithm(Parameter &param)
{
    mpf_class r;
    mpf_class a = (param.a + param.b) / 2;
    mpf_class b = sqrt(param.a * param.b);
    mpf_class t = param.t - param.p * power((param.a - a), 2);
    mpf_class p = 2 * param.p;

    r = power((a + b), 2) / (4 * t);

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

    return r;
}

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

    //★精度 ここを調整する
    int digits = 200;
    //
    //精度
    mpf_set_default_prec(digits * log2(10));
    //何回計算する？
    int recursion_count = log2(digits);

    Parameter param;
    mpf_class pi;

    for (int i = 1; i <= recursion_count; i++)
    {
        pi = gauss_Legendre_algorithm(param);
    }

    std::cout << "有効桁 = " << digits << std::endl;
    std::cout << "計算回数 = " << recursion_count << std::endl;

    mp_exp_t exp;
    std::string significand = pi.get_str(exp);
    std::cout << "pi = " << significand.insert(1, ".") << std::endl;

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

結果

Start..
有効桁 = 200
計算回数 = 7
pi = 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097565
End..

計算結果
3.
1415926535 8979323846 2643383279 5028841971 6939937510
5820974944 5923078164 0628620899 8628034825 3421170679
8214808651 3282306647 0938446095 505822

実際は
3.
1415926535 8979323846 2643383279 5028841971 6939937510
5820974944 5923078164 0628620899 8628034825 3421170679
8214808651 3282306647 0938446095 5058223172 5359408128

おお、当たり前かもしれないけど合っている！！。

ちなみに　int digits = 100万(1000000)で計算してみる。
計算回数はたったの19回、計算時間は数秒！！

結果

Start..
有効桁 =1000000
計算回数 = 19
pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865
...長いので省略！！
403332127284919441843715069655208754245059895678796130331164628399634646042209010610577945815130927562832084527
End..