Help us understand the problem. What is going on with this article?

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

More than 3 years have 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時間もバグ探しをしてしまいました。(汗

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away