円周率646456700桁弱を計算する #1

gmpの限界を見ておこうと思う。環境は以下の通り。
gmp-5.0.5 /gcc-4.4.5 /debian64-2.6.32-5-amd64 /VirtualBox 4.2.4 /Windows 7 /DELL XPS8300 /3.4G Intel Core i7-2600。

「GNU/Linux上で円周率の計算をおこなう 」すずきひろのぶ氏の百万桁計算のソースを使用する。変更点は以下の通り。itchynyさんの修正も適用している。すずきひろのぶ氏が指摘している通り、gmp の precision のサイズの上限により、円周率計算は6.4億桁が上限となる。12/15追記：大嘘です。円周率10億桁を計算する #1

nitobe@debian64:~/pi$ diff pi.c pi2.c
16a17
>   mpz_t temp, temp2;/* add itchyny */
18c19
<   dprec = 1000000L;/* decimal precision */
---
>   dprec = 800000000L;/* decimal precision */
21c22
<   loopmax = 21;
---
>   loopmax = 31;
22a24
>   fprintf(stderr, "%d initial:", loopmax); /* Add nitobe */
34,35c36,38
< 
<    mpf_set_ui (A, 1);
---
>   mpz_init (temp);/* add itchyny */
>   mpz_init (temp2);/* add itchyny */
>   mpf_set_ui (A, 1);
41a45
>   mpz_set_ui (temp2, 2);/* add itchyny */
43a48
>       fprintf(stderr, " %d:", k); /* Add nitobe */
53c58,60
<       mpf_pow_ui (t2, c2, k);/* 2^k */
---
> //      mpf_pow_ui (t2, c2, k);/* 2^k */
>       mpz_pow_ui (temp, temp2, k);/* 2^k */
>       mpf_set_z (t2, temp);/* add itchyny */
59c66,67
<   mpf_pow_ui (t2, t1, 2);
---
> //  mpf_pow_ui (t2, t1, 2);
>   mpf_mul (t2, t1, t1);/* mod itchyny */
64c72
<   exit(0);
---
>   return 0;/* mod nitobe */
nitobe@debian64:~/pi$ 

dprecで求める桁数を指定している。 loopmaxは演算繰り返し回数だが、求め得る桁数と、以下の関係がある。

loopmax	correct digits (about)
9	<500
10	<800
11	<1500
12	<2900
13	<5700
14	<11300
15	<22500
16	<44900
17	<89600
18	<179000
19	<357800
20	<715500
21	<1430800
22	<2861500
23	<5722800
24	<11445400
25	<22890600
26	<45781000
27	<91561900
28	<183123600
29	<366247100
30	<732494000

makefile

nitobe@debian64:~/pi$ cat makefile
CC=gcc
CFLAG=-static -O3
GMP=gmp-5.0.5
GMP_INC=-I/usr/local/$(GMP)/include
GMP_LIB=/usr/local/$(GMP)/lib/libgmp.a

.PHONY: clean

all: pi

pi:
        $(CC) $(CFLAG) $(GMP_INC) pi2.c -o pi $(GMP_LIB)

clean:
rm pi
nitobe@debian64:~/pi$

make して実行

nitobe@debian64:~/pi$ make
gcc -static -O3 -I/usr/local/gmp-5.0.5/include pi2.c -o pi /usr/local/gmp-5.0.5/lib/libgmp.a
nitobe@debian64:~/pi$ time ./pi > pi800M.txt
31 initial: 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19: 20: 21: 22: 23: 24: 25: 26: 27: 28: 29: 30: 31:
real68m26.877s
user66m50.763s
sys1m36.338s
nitobe@debian64:~/pi$

単一行巨大ファイルは扱い難いので、百桁／行に整形する。
format.sed

nitobe@debian64:~/pi$ cat format.sed
#!/bin/sed -f
/^$/d
s/0.3/3.¥n/
s/[0-9]¥{100¥}/&¥n/g
s/e1/¥n/
nitobe@debian64:~/pi$ time ./format.sed pi800M.txt >pi800m.txt

real1m1.508s
user1m0.404s
sys0m1.108s
nitobe@debian64:~/pi$ ls -al
合計 1271884
drwxr-xr-x  2 nitobe nitobe      4096 2012-12-13 21:15 .
drwxr-xr-x 37 nitobe nitobe      4096 2012-12-13 19:50 ..
-rwxr-xr-x  1 nitobe nitobe        62 2012-12-04 21:16 format.sed
-rw-r--r--  1 nitobe nitobe       208 2012-12-13 19:49 makefile
-rwxr-xr-x  1 nitobe nitobe    854298 2012-12-13 19:49 pi
-rw-r--r--  1 nitobe nitobe      1851 2012-12-13 17:05 pi.c
-rwxr-xr-x  1 nitobe nitobe    854298 2012-12-13 19:40 pi2
-rw-r--r--  1 nitobe nitobe      2264 2012-12-13 19:48 pi2.c
-rw-r--r--  1 nitobe nitobe 646456999 2012-12-13 20:58 pi800M.txt
-rw-r--r--  1 nitobe nitobe 652921567 2012-12-13 21:08 pi800m.txt
nitobe@debian64:~/pi

８億桁を指定しているが、６億４千万桁しかない。フォーマット変換後、改行「¥n」が６百万個増えている。百桁／行フォーマットなのでエディタで開ける。重いけど。これは vi で見てみた。

      1 3.
      2 1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679
      3 8214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196
　　　　：
　　　　：
6464569 9336326991629433891751838174109742989418332960270103331518724202832616454999878479262174537890538106
6464570 6455256718445892438033533034196466727204878408211469825653223223178441813002292106453186861390731397
6464571 2309985607826687727041129258325166030790611057861923604585439381176655074474471470989090334958

さて、この値が正しいかどうか？それが問題だ。つづく

pi