円周率πを100万桁まで計算するプログラムを作ってみました。super pi もどきです。ソースコードと解説を載せます。
ソースコードは80行程度の非常に短いものですが、このような簡単な実装でも最近のPCなら100万桁の円周率が一分も掛らず計算できます。
ソースコード
/** * pi.cpp * * written by pyopyopyo at gmail dot com */ #include <math.h> #include <gmpxx.h> #include <iostream> // N桁の円周率を計算 #define N 100*100*100 using namespace std; int main() { // N桁の数値計算に必要な精度を prec(bits) とすると // // prec = log2(10^N) // = N * log2(10) // 計算誤差を考えるのがめんどうなので，とりあえず(N+5)桁で計算してみる int prec = (N+5)*log2(10); mpf_set_default_prec(prec); fprintf(stderr, "prec:= %d bits\n", prec); const mpf_class two ("2.0"); mpf_class a ("1.0"); mpf_class b (a/sqrt(mpf_class("2.0"))); mpf_class t (a); mpf_class x (a); mpf_div_ui(t.get_mpf_t(), t.get_mpf_t(), 4); mpf_class tmp(a); const int P = ceil(log2(N)); for (int k=1; k<=P; ++k) { tmp = a; a = a + b; mpf_div_ui(a.get_mpf_t(), a.get_mpf_t(), 2); b = sqrt( tmp*b ); tmp -= a; mpf_pow_ui(tmp.get_mpf_t(), tmp.get_mpf_t(), 2); t -= x * tmp; x *= two; } tmp = a + b; mpf_pow_ui(tmp.get_mpf_t(), tmp.get_mpf_t(), 2); tmp /= t; mpf_div_ui(tmp.get_mpf_t(), tmp.get_mpf_t(), 4); 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, &str[i], 10); cout << tmp ; if (0 == (count%10)) cout << endl; else cout <<" "; } return 0; }
コンパイル方法と、実行方法
$ gcc -lgmpxx pi.cpp -o pi
$ time ./pi
計算時間の目安としては
ぐらいになります。
なお出力は、以下のURLで公開されている１００万桁の円周率πのデータの形式に合わせています。
解説
- 計算精度の指定
int prec = (N+5)*log2(10); mpf_set_default_prec(prec);
gmpのデフォルトの計算精度を、(N+5)桁に設定しています。
log2(10^(N+5)) → (N+5)*log2(10) ということです。
なお、今回はN=100万なので、 precは約332万になります。つまり、GMPの変数は、一変数につき 400Kbyte 程度のメモリを消費することになります。
gmpのc++インタフェースを利用します。ヘッダファイルは gmpxx.h です。
#include <gmpxx.h>
c++インタフェースを使うと、たとえば、gmpのcインタフェースで書いた以下のコードは
mpf_t a, b, c; mpf_init(a); mpf_init(b); mpf_init(c); mpf_set_str(a, "1234", 10); mpf_set_str(b, "5678", 10); mpf_add(c, a, b); mpf_out_str(stdout, 10, 10, c);
c++インタフェースでは
mpf_class a,b,c; a="1234"; b="5678"; c=a+b; cout << c;
と簡潔に記述できます。
さらに
mpf_class a("1234"),b("5678"),c; c=a+b; cout << c;
とも書け、さらに圧縮して
cout << mpf_class("1234")+mpf_class("5678");
と書くことも出来ます。
おすすめ度の平均: 
意外と忘れがちだったり、重要だけど知らなかったりするのが 常識なんです。 文字どおり「常識」にしたい内容 困る前に
ニューメリカルレシピ・イン・シー 日本語版―C言語による数値計算のレシピ
posted with amazlet at 09.03.08
