LoginSignup
2

More than 1 year has passed since last update.

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

Last updated at Posted at 2017-12-11

初めに

今回は任意精度数演算ライブラリ、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  #エラー

結果

一応、円周率が適切に求められました。
結果が長すぎるので流石に載せません。

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
2