Posted at

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

More than 1 year has passed since last update.


初めに

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


結果

一応、円周率が適切に求められました。

結果が長すぎるので流石に載せません。


ちょっとした悲劇

今回円周率を求めるプログラムを実装するにあたって、初心者にありそうな悲劇が起こりました。

上記のプログラムなのですが最初は、19行目を

mpf_class t(1/4,65546);

と書いてしまったのです。

これだと、C言語では整数割る整数なので自動的に切り捨ててしまい、0になってしまいます。

これを見つけるのに1時間もバグ探しをしてしまいました。(汗