多倍長計算ライブラリGMPで円周率を求める。
多倍長計算したいなー、と思ってググってみるとGMPなるライブラリがあったので、試してみる。
触り方も何も分からんかったので、この辺を参照。
どうやら
mpf_t 多倍長変数(構造体?)を宣言し
mpf_init で変数の初期化をして
以下の関数で
mpf_add 第一引数に第二引数と第三引数の和を
mpf_sub 同様に差を
mpf_mul 同様に積を
mpf_div 同様に商を代入する
らしい。
どうせなので参考ページとは変えてガウス=ルジャンドルのアルゴリズムで組んでみる。
#include#include #define LOG_TEN_TWO 3.32192809488736234789 #define bprec(n) (int)(((n+10)*LOG_TEN_TWO)+2) int main(int argc,char *argv[]) { mpf_t a,b,t,x,y,p,t1,t2,t3,t4; long i; long int prec, dprec; dprec = 1000000L; prec = bprec(dprec); mpf_set_default_prec (prec); mpf_init(a); mpf_init(b); mpf_init(t); mpf_init(x); mpf_init(y); mpf_init(p); mpf_init(t1); mpf_init(t2); mpf_init(t3); mpf_init(t4); mpf_set_ui(a,1); //a = 1 mpf_set_ui(t1,1); mpf_div_ui(t2,t1,2); mpf_sqrt(b,t2); //b = sqrt(1/2) mpf_div_ui(t,t1,4); //t = 1/4 mpf_set_ui(p,1); //p = 1 for(i = 0;i<10;i++) { mpf_add(t1,a,b); mpf_div_ui(x,t1,2); //x = (a+b)/2 mpf_mul(t1,a,b); mpf_sqrt(y,t1); //y = sqrt(a*b) mpf_sub(t1,a,x); mpf_pow_ui(t2,t1,2); mpf_mul(t3,p,t2); mpf_sub(t4,t,t3); mpf_set(t,t4); //t = t-p(a-x)^2 mpf_set(a,x); //a = x mpf_set(b,y); //b = y mpf_mul_ui(t1,p,2); mpf_set(p,t1); //p = 2*p } mpf_add(t1,a,b); mpf_mul(t2,t1,t1); mpf_mul_ui(t3,t,4); mpf_div(t4,t2,t3); mpf_out_str(stdout,10,1000000,t4); printf("\n"); return 0; }
調子に乗って100万桁とかにしているが、ループ回数が少なすぎて、そこまで精度が出ていない気がするorz
(出力はしっかり100万桁出しているが。あとでちゃんとトレーススクリプト書いてみようっと。・・・2786桁しかあってなかったorzもっとループを増やせば!)
とりあえずGMPの概要が分かったので、ここまで。
CC = gcc GMP = /usr/lib/libgmp.a pi:pi.c $(GMP) $(CC) -o $@ $< $(GMP)