円周率を求めるプログラムを書く

初めに

今回は任意精度数演算ライブラリ、GMP(GNU Multi-Precision Library)を使って円周率を求めたいと思います。
円周率を求めるアルゴリズムはガウス＝ルジャンドルのアルゴリズムを使います。

ソースコード

main.cpp

#include <cstdio>
#include <iostream>
#include <gmp.h>
#include <gmpxx.h>
#include <cmath>
//#include <omp.h>


int main(){

    int n=0;
    std::cout<<"n?>";
    std::cin>>n;

    mpf_class a(1,65546);//初期値1,65546bit
    mpf_class ba(0,65546);
    mpf_class b(1.0/sqrt(mpf_class(2)),65546);
    mpf_class t(1.0/4.0,65546);
    mpf_class p(1,65546);

    //openMP
    //#pragma omp parallel for
    for(int i=0;i<n;i++){
        ba=a;
        a=(ba+b)/2;
        b=sqrt(ba*b);
        t=t-(p*(ba-a)*(ba-a));
        p=2*p;
    };

    mpf_class A=(a+b)*(a+b)/(4*t);

    mp_exp_t exp;

    std::string A_str=A.get_str(exp);

    std::cout<<A_str.insert(1, ".")<<std::endl;
}

アルゴリズムはwikipediaを見ながら実装しました。
変数のビット長は65546bitにしています。
速度向上のため、openMPを使おうとしたのですが漸化式の形になっているので上手く速度が上がりませんでした。

コンパイル

g++ main.cpp -lgmp -lgmpxx

下記のリンクの順番だと私の環境(OS、Windows 10。コンパイラ、mingw)では上手くいきませんでした。

g++ -lgmp -lgmpxx main.cpp