( 続 ) 円周率を求める方法

前回のプログラムを C++ で書きなおして任意精度演算を行ってみました。

結果

だいたい1000回ループを回して小数点以下599桁まであてました。600桁で2のところが1になってました。惜しかったです。

感想

自分で考えたアルゴリズムがきちんと動作してて数学の確かさを確認するのと同時に感動しました。円周率がきちんと自分で求められてしかも無限回すれば正確な値が出るのはなかなか達成感があります(無限にできないしpiは無限小数だから真に正確な値はでませんが)。

コード

一応使ったコード置いときます。C++はほとんど触ったことがないのでsuper pi を作ろうを参考にパクって自分のコードが動くよう直しました。

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

using namespace std;

#define  N  10000

int main() {
    int prec = (N+1)*log2(10);
    mpf_set_default_prec(prec);

    mpf_class a ("1.0");
    mpf_class b (a/sqrt(mpf_class("2.0")));
    mpf_class c (b);
    mpf_class _a,_b;
    int m = 1000;
    for(int n = 0; n < m; n++) {
        _b = b*(c/(a+c));
        _a = sqrt(0.5+_b*_b);
        b = _b;
        a = _a;アルゴリズム
    }
    mpf_class tmp (b/a*pow(2,m-1)*8);

    //すいません! 以下コピペです
    mp_exp_t expo;
    char *str = mpf_get_str(0, &expo, 10, N+1, tmp.get_mpf_t());

    for (int i=expo,count=1; i<N; i+=10, ++count) {
        char tmp[12]={0};
        strncpy(tmp, &amp;str[i], 10);
        cout << tmp ;
        if (0 == (count%10))
            cout << endl;
        else
            cout << " ";
    }
    //ここまでコピペです。すいません。
    cout << endl;
    return 0;
}