円周率10億桁を計算する #1

---

すずきひろのぶ氏が指摘している通り、gmp の precision のサイズの上限により、
円周率計算は6.4億桁が上限となる。

m(_ _)m　大変申し訳ございません。嘘つきました。勉強不足。
現在、gmp の precision のサイズは int (2^31-1=2,147,483,647)ではありません。従って、(2^31-1)/log210 ≑ 646,456,992 が円周率計算の上限ではありません。

nitobe@debian64:~/pi$ cat size.c
#include<stdio.h>
#include<gmp.h>

main()
{
printf("int: %d¥n", sizeof(int));
printf("long int: %d¥n", sizeof(long int));
printf("long long int: %d¥n", sizeof(long long int));
printf("mp_bitcnt_t: %d¥n", sizeof(mp_bitcnt_t));
}
nitobe@debian64:~/pi$ make size
gcc -static -O3 -m64 -mtune=amdfam10 -I/usr/local/gmp-5.0.5/include size.c -o size /usr/local/gmp-5.0.5/lib/libgmp.a
nitobe@debian64:~/pi$ ./size
int: 4
long int: 8
long long int: 8
mp_bitcnt_t: 8
nitobe@debian64:~/pi$ 

理論上の上限は、2^63-1 = 2進 9,223,372,036,854,775,807桁、(2^63-1)/log210 = 10進 2,776,511,644,261,678,565桁 ということになる。2.78E(エクサ10^18)Bit = 348P(ペタ10^15)Byteのメモリを実装したコンピュータなどそうそうあるもんじゃない。当面 gmp の上限は見られないということだ。

鋭意計算中。

12/15 23:40 追記：
pi2.c　(diff)

nitobe@debian64:~/pi$ diff pi.c pi2.c
10,11c10,11
< #define LOG_TEN_TWO  3.32192809488736234789
< #define bprec(n) (int)(((n+10)*LOG_TEN_TWO)+2)
---
> #define LOG_TEN_TWO  3.3219280948873623478703194294894
> #define bprec(n) (long int)(((n+10)*LOG_TEN_TWO)+2)
16a17
>   mpz_t temp, temp2;/* add itchyny */
18c19
<   dprec = 1000000L;/* decimal precision */
---
>   dprec = 1000000000L;/* decimal precision */
21c22
<   loopmax = 21;
---
>   loopmax = 31;
22a24
>   fprintf(stderr, "%ld %d initial:", prec, 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);
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 */
61a70
>   fprintf(stderr, " conversion:", prec, loopmax); /* Add nitobe */
64c73
<   exit(0);
---
>   return 0;/* mod nitobe */
nitobe@debian64:~/pi$ 

makefile

nitobe@debian64:~/pi$ cat makefile
CC=gcc
OPT= -static -O3 -m64 -mtune=amdfam10
GMP=gmp-5.0.5
GMP_INC=-I/usr/local/$(GMP)/include
GMP_LIB=/usr/local/$(GMP)/lib/libgmp.a
SRC= pi2.c
EXE= pi

.PHONY: clean

all: $(EXE)

pi: $(SRC)
        $(CC) $(OPT) $(GMP_INC) $(SRC) -o $(EXE) $(GMP_LIB)

size: size.c
        $(CC) $(OPT) $(GMP_INC) size.c -o size $(GMP_LIB)

clean:
rm -f pi
rm -f size
nitobe@debian64:~/pi$ 

nitobe@debian64:~/pi$ time ./pi >pi1G.txt
3321928130 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: conversion:
real    196m38.626s
user    192m55.963s
sys     3m43.246s
nitobe@debian64:~/pi$ time ./format.sed pi1G.txt >pi1g.txt

real2m54.283s
user2m39.994s
sys0m3.148s
nitobe@debian64:~/pi$

注：実行時間は全く当てにならない。このPC上では、seti@home on BOINC 及び、World Community Grid on BOINC がリソースの75%を占有している。

pi1g.txt

       1 3.
       2 1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679
       3 8214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196

          :
          :
10000000 4219776753871319682188195635848934815504410194647387557034502943416861599324354199731814355060392734
10000001 6434543524276655356743570219396394581990548327874671398682093196353628204612755715171395115275045519
10000002 64388312

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

pi