理系の憧れ！円周率の求め方～long double の壁を乗り越える～

はじめに

この記事は前回のプログラムを改良したものについて扱います。ソースコードの意味を知りたい方は、前回の記事を先にご覧ください。

前回までの方法では、いくら正確に円周率を計算しても、long doubleの精度が限界でした。そこで、精度を理論上、メモリの限界まで上げることができるライブラリである、"GMP"というものを使用しました。これにより、8GBメモリ環境では余裕で前回の精度を超えることができました。

GMPとは

GMPは"GNU Multiple Precision library"の略で、多倍長の計算を簡単に実装できるライブラリとなっています。ちなみに、これはUnix系のOSでのみ動作します。Windowsで動かすには、Cygwinを入れるか、GMPの代わりにMPIRというものを入れる必要があります。幸運にも、私のPCは、LinuxとWindowsをデュアルブートしていたので、Linux環境で計算しました。Windows環境をご利用の方は、Cygwinを入れるのをお勧めします。(Linuxのデュアルブート環境を構築するのもアリですが、いろいろ面倒なことになるので、上級者向けです…)

環境

OS:Ubuntu 16.04

Editor:Visual Studio Code

CPU:Intel Core i7-2600K 3.40GHz

RAM:8.00GB

ちなみに、これはサブPCです。。。

GMPに移植

この作業が本当に大変でした…。

初期化する関数を使わなければならない。
四則演算はすべて専用の関数で行う。+-*/は使えない。
メモリを最後に開放しなければならない。
計算ごとに中間状態を保存するtemp変数を使わなければならない。

最後、どうしてもセグフォになると思って、ググってみると、変数の初期化、開放関数を複数の変数に対して同時に使う時は、最後の引数にNULLを入れないといけないとか。~~そんなの知らんぞ。。。~~自分の調べた限りでは、英語のサイトしか出てこなかったので、余計厄介でした。ちゃんとリファレンスを読めってことですね。

コードを移植した結果

こうなりました。


#include <math.h>
#include <gmp.h>
#include <stdio.h>

int main(void) {
	mpz_t n, i, temp1z;
	mpf_t a,b,t,temp1,temp2,temp4,p;
	mpf_set_default_prec(10000);
	mpz_init_set_str(n,"100",10);//回数
	mpz_init_set_str(temp1z,"1",10);
	mpf_init_set_str(p,"1",10);
	mpz_init_set_str(i,"0",10);
	mpf_init_set_str(a,"1",10);
	mpz_init_set_str(i,"0",10);
	mpf_init_set_str(temp1,"1",10);
	mpf_init_set_str(temp2,"2",10);
	mpf_init_set_str(temp4,"4",10);
	mpf_init(b);
	mpf_init(t);

	mpf_sqrt(b,temp2);
	mpf_div(b,b,temp2);
	mpf_div(t,a,temp4);
	
	while(mpz_cmp(n,i)){
		mpf_t a2,b2,t2,tempx,tempy,tempz,tempa,tempb;
		mpf_init(a2);
		mpf_init(b2);
		mpf_init(t2);
		mpf_init(tempx);
		mpf_init(tempy);
		mpf_init(tempz);
		mpf_init(tempa);
		mpf_init(tempb);
		mpf_add(tempx,a,b);
		mpf_div(a2,tempx,temp2);
		mpf_mul(tempy,a,b);
		mpf_sqrt(b2,tempy);
		mpf_sub(tempz,a,a2);
		mpf_mul(tempa,tempz,tempz);
		mpf_mul(tempb,p,tempa);
		mpf_sub(t2,t,tempb);
		mpf_mul(p,p,temp2);
		mpf_set(a,a2);
		mpf_set(b,b2);
		mpf_set(t,t2);
		mpz_add(i,i,temp1z);
		mpf_clear(a2);
		mpf_clear(b2);
		mpf_clear(tempb);
		mpf_clear(tempa);
		mpf_clear(tempz);
		mpf_clear(tempy);
		mpf_clear(tempx);
		mpf_clear(t2);
	}

	mpf_t tempu,tempd,ans;
	mpf_inits(tempu,tempd,ans,NULL);
	mpf_add(tempu,a,b);
	mpf_mul(tempu,tempu,tempu);
	mpf_mul(tempd,temp4,t);
	mpf_div(ans,tempu,tempd);

	FILE *text;
	text = fopen("pi","w+");
	mpf_out_str(text,10,1000,ans);
	fclose(text);
	
	mpz_clears(n,i,NULL);
	mpf_clears(p,a,b,t,temp1,temp2,temp4,tempu,tempd,ans,NULL);

	return 0;
}

やっていることは、前回の記事とほぼ同じです。複雑に見えますが。p2を経由しないように変更したり、出力先をテキストファイルにしたり。GMPに興味のある方は、こんなサイトなどを参考にしてみてください。

出力

0.3141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420199e1

これが生のデータです。最後についているe1は、10の1乗という意味なので、正しく3.14...になっていることがわかります。さて、天下のWikiから、正しいπのデータを引用してきます。

以下1000桁までの値
$π$ = 3.
14159 26535 89793 23846 26433 83279 50288 41971 69399 37510 58209 74944 59230 78164 06286 20899 86280 34825 34211 70679 82148 08651 32823 06647 09384 46095 50582 23172 53594 08128 48111 74502 84102 70193 85211 05559 64462 29489 54930 38196 44288 10975 66593 34461 28475 64823 37867 83165 27120 19091 45648 56692 34603 48610 45432 66482 13393 60726 02491 41273 72458 70066 06315 58817 48815 20920 96282 92540 91715 36436 78925 90360 01133 05305 48820 46652 13841 46951 94151 16094 33057 27036 57595 91953 09218 61173 81932 61179 31051 18548 07446 23799 62749 56735 18857 52724 89122 79381 83011 94912 98336 73362 44065 66430 86021 39494 63952 24737 19070 21798 60943 70277 05392 17176 29317 67523 84674 81846 76694 05132 00056 81271 45263 56082 77857 71342 75778 96091 73637 17872 14684 40901 22495 34301 46549 58537 10507 92279 68925 89235 42019 95611 21290 21960 86403 44181 59813 62977 47713 09960 51870 72113 49999 99837 29780 49951 05973 17328 16096 31859 50244 59455 34690 83026 42522 30825 33446 85035 26193 11881 71010 00313 78387 52886 58753 32083 81420 61717 76691 47303 59825 34904 28755 46873 11595 62863 88235 37875 93751 95778 18577 80532 17122 68066 13001 92787 66111 95909 21642 01989 …

比較してみると、正しく1000桁(正確には998桁)計算できていることがわかります。

プログラムのパラメータを書き換えれば、メモリの限り精度をいくらでも上げられます。ちなみに、1000桁程度では、enterを叩いた瞬間に計算が終わるほどの速さです。

また機会があれば、高みを目指していきたいものです。ご覧いただき、ありがとうございました。

この記事への感想を教えてください

内容が十分 ()
内容が足りなかったが役立った ()
内容が足りず役立たなかった ()
求めている記事ではなかった ()