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


2012-12-18 付で、多倍長演算ライブラリ gmpLink が、5.1.0 になった。

早速buildしてみよう。環境は、gcc-4.4.5 /debian64-2.6.32-5-amd64 /VirtualBox 4.2.4 /Windows 7 /DELL XPS8300 /3.4G Intel Core i7-2600。ramは、16GByte実装、うち12GByteをVirtualBoxに割り当てている。
root@debian64:/home/nitobe# wget ftp://ftp.gmplib.org/pub/gmp-5.1.0/gmp-5.1.0.tar.bz2Link 
     :
root@debian64:/home/nitobe# tar xvf gmp-5.1.0.tar.bz2
     :
root@debian64:/home/nitobe# cd gmp-5.1.0
root@debian64:/home/nitobe# mkdir build
root@debian64:/home/nitobe# cd build
root@debian64:/home/nitobe# ../configure -prefix=/usr/local/gmp-5.1.0
     :
root@debian64:/home/nitobe# make -j 4
     :
root@debian64:/home/nitobe# make check
     :
root@debian64:/home/nitobe# make install
     :
root@debian64:/home/nitobe# tree -L 1 /usr/local
/usr/local
├── bin
├── etc
├── games
├── gmp-4.2.2
├── gmp-5.0.5
├── gmp-5.1.0
├── include
├── lib
├── man -> share/man
├── sbin
├── share
└── src

12 directories, 0 files
root@debian64:/home/nitobe#

makefile を整えて
nitobe@debian64:/home/nitobe/itchyny# cat makefile
CC=gcc
OPT= -static -O3 -m64 -mtune=amdfam10
GMP422=gmp-4.2.2
GMP422_INC=-I/usr/local/$(GMP422)/include
GMP422_LIB=/usr/local/$(GMP422)/lib/libgmp.a
GMP505=gmp-5.0.5
GMP505_INC=-I/usr/local/$(GMP505)/include
GMP505_LIB=/usr/local/$(GMP505)/lib/libgmp.a
GMP510=gmp-5.1.0
GMP510_INC=-I/usr/local/$(GMP510)/include
GMP510_LIB=/usr/local/$(GMP510)/lib/libgmp.a
SRC= pi.c
EXE= pi

.PHONY: clean

all: $(EXE)

pi: $(SRC)
        $(CC) $(OPT) $(GMP510_INC) $(SRC) -o $(EXE) $(GMP510_LIB)

pi422: $(SRC)
        $(CC) $(OPT) $(GMP422_INC) $(SRC) -o $(EXE)422 $(GMP422_LIB)

pi505: $(SRC)
        $(CC) $(OPT) $(GMP505_INC) $(SRC) -o $(EXE)505 $(GMP505_LIB)

pi510: $(SRC)
        $(CC) $(OPT) $(GMP510_INC) $(SRC) -o $(EXE)510 $(GMP510_LIB)

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

clean:
        rm -f pi
        rm -f pi422
        rm -f pi505
        rm -f pi510
        rm -f size
nitobe@debian64:/home/nitobe/itchyny# 
まずは、円周率1億桁の実行時間を比較する。使用するソースは itchyny Link さんのもの。itchyny さんのソースをちょっと改変。itchyny さん、ごめんなさい。
nitobe@debian64:/home/nitobe/itchyny# diff pi.c pi.c.org
66,72c66,68
<   long int digits = 100000000;
<   long int prec = digits * log2(10);
<   long int n = digits / 14;
< 
<   fprintf(stderr, "%ld digits(prec=%ld) of pi, %ld terms.¥n", digits, prec, n);
<   fprintf(stderr, "with gcc-%d.%d.%d, gmp-%s¥n",
<                 __GNUC__, __GNUC_MINOR__, __GNUC_PATCHLEVEL__, gmp_version);
---
>   int digits = 10000;
>   int prec = digits * log2(10);
>   int n = digits / 14;
125c121
<   mpf_out_str(fp, 10, digits + 10, pi);
---
>   mpf_out_str(fp, 10, digits, pi);
142d137
< 
nitobe@debian64:/home/nitobe/itchyny#
それぞれの gmp バージョンの実行ファイルをコンパイルして、実行。
nitobe@debian64:/home/nitobe/itchyny# make pi422
gcc -static -O3 -m64 -mtune=amdfam10 -I/usr/local/gmp-4.2.2/include pi.c -o pi422 /usr/local/gmp-4.2.2/lib/libgmp.a
nitobe@debian64:/home/nitobe/itchyny# make pi505
gcc -static -O3 -m64 -mtune=amdfam10 -I/usr/local/gmp-5.0.5/include pi.c -o pi505 /usr/local/gmp-5.0.5/lib/libgmp.a
nitobe@debian64:/home/nitobe/itchyny# make pi510
gcc -static -O3 -m64 -mtune=amdfam10 -I/usr/local/gmp-5.1.0/include pi.c -o pi510 /usr/local/gmp-5.1.0/lib/libgmp.a
nitobe@debian64:/home/nitobe/itchyny$ time ./pi422
with gcc-4.4.5, gmp-4.2.2
digits=100000000; prec=332192809; n=7142857
356.010 s
491.560 s

real8m11.823s
user8m9.895s
sys0m1.692s
nitobe@debian64:/home/nitobe/itchyny$ time ./pi505
with gcc-4.4.5, gmp-5.0.5
digits=100000000; prec=332192809; n=7142857
186.400 s
230.130 s

real3m50.405s
user3m48.514s
sys0m1.648s
nitobe@debian64:/home/nitobe/itchyny$ time ./pi510
with gcc-4.4.5, gmp-5.1.0
digits=100000000; prec=332192809; n=7142857
162.630 s
200.160 s

real3m20.511s
user3m18.580s
sys0m1.616s
nitobe@debian64:~/itchyny$
4.2.2~5.0.5:63%、5.0.5~5.1.0:14% の速度向上だ。
実は、上の改変で、10億桁も計算可能になっている。66行目を 1000000000 にして、make
nitobe@debian64:~/itchyny$ make pi510
gcc -static -O3 -m64 -mtune=amdfam10 -I/usr/local/gmp-5.1.0/include pi.c -o pi510 /usr/local/gmp-5.1.0/lib/libgmp.a
nitobe@debian64:~/itchyny$ time ./pi510
1000000000 digits(prec=3321928094) of pi, 71428571 terms.
with gcc-4.4.5, gmp-5.1.0
2501.230 s
3059.550 s

real51m26.963s
user50m18.689s
sys0m41.095s
nitobe@debian64:~/itchyny$ cat format.sed
#!/bin/sed -f
/^$/d
s/^0.3/3.¥n/
s/[0-9]¥{100¥}/&¥n/g
s/e1//
nitobe@debian64:~/itchyny$ time ./format.sed output.txt >pi1g.txt

real2m4.889s
user1m38.666s
sys0m1.520s
nitobe@debian64:~/itchyny$ diff pi1g.txt ../pi/pi1g.txt
10000002c10000002
< 64388312
¥ ファイル末尾に改行がありません
---
> 64388312
nitobe@debian64:~/itchyny$ 

円周率100万桁のベンチマークは、「円周率ベンチマークをやり直すLink 」を参照していただきたい。5.0.5と比較して、35%~25%程度速くなっている。何故か、2.0G Intel(R) Core(TM) Duo CPU T2500 と Intel(R) Pentium(R) M 1.4GHz は速度の向上が全く見られない。1.2 GHz ARM Marvell Kirkwood 88F6281 SheevaPlug は、makeでエラーとなりbuild 不能。

現在の環境で、頑張れば20億桁位まで行けそうだけど、あとは ram の実装次第だよね。64GByte実装で100億桁は行けるかな?さもなくば外部記憶に逃がすか?遅くなるなぁ。

d72a1888159a65c16f1b084b477fa18b

Moore's Law (集積回路上のトランジスタ数は「18か月ごとに倍になる」)を信頼して、暫し待つか。

 


— posted by nitobe at 11:11 am   commentComment [1] 

ちはやふる 19



51M6Osnr9ML__SS500_


ちはやふる(19) (BE LOVE KC)
著者: 末次 由紀
出版社: 講談社
ISBN: 4063803694
定 価: ¥ 450
出版日: 2012-12-13

本宅には使えなかった。

 

— posted by nitobe at 09:41 pm   commentComment [0] 

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


と、いうわけで、検証。
macroxue / const-piLink のzipファイルを入手する。
件の Hanhong Xue 氏の円周率計算プログラム、新しめのバージョンである。
nitobe@debian64:~$ wget https://github.com/macroxue/const-pi/archive/master.zipLink 
     :
     :
nitobe@debian64:~$ unzip master.zip
Archive:  master.zip
8232a2005009a1d3f33c68c2e7e1054d12e8ae6d
   creating: const-pi-master/
  inflating: const-pi-master/Args.h  
  inflating: const-pi-master/README.md  
  inflating: const-pi-master/const.cpp  
  inflating: const-pi-master/fac.c   
  inflating: const-pi-master/fac.h   
  inflating: const-pi-master/makefile  
  inflating: const-pi-master/my.c    
  inflating: const-pi-master/my.h    
nitobe@debian64:~$ 
10億桁の計算はできるのだが、2進10進変換で落ちる。たぶんメモリ不足だ。フォーマットも非常に凝ったものだが、扱いにくいので、gmpのオリジナル出力関数を使用する。my.c の my_out_str関数もコメントアウトして亡き者とする。Xue さんごめんなさい。
nitobe@debian64:~/const-pi-master$ diff const.cpp ../const_new/const-pi-master/const.cpp
90c90,91
<             my_out_str(out_fp, 10, digits+2, result);
---
> //            my_out_str(out_fp, 10, digits+2, result);
>             mpf_out_str(out_fp, 10, digits+2+2, result);
nitobe@debian64:~/const-pi-master$ diff my.c ../const_new/const-pi-master/my.c
367a368
> /*
434a436
> */
ささ、早速計算してみよう。
nitobe@debian64:~/const_new/const-pi-master$ time ./const pi -d1000000000 -o pi1G.txt
Calculation started at Sat Dec 15 12:47:03 2012

1000000000 digits of pi using Chudnovsky formula, 70513670 terms.
        initialization : 17.720 s
      binary splitting : 4126.380 s                
              division : 771.880 s
           square root : 115.470 s
        multiplication : 76.110 s
-----------------------------------
         calculation : 5108.710 s
              output : 2377.740 s

           user time : 7453.158 s
            sys time : 33.330 s
        memory usage : 7567424 KB

Calculation finished at Sat Dec 15 14:51:49 2012

real124m46.183s
user124m13.158s
sys0m33.342s
nitobe@debian64:~/const_new/const-pi-master$ time ./format.sed pi1G.txt >pi1g.txt

real2m36.393s
user2m34.690s
sys0m1.696s
nitobe@debian64:~/const_new/const-pi-master$ time ./const pi-rama -d1000000000 -o pi1G-rama.txt
Calculation started at Sat Dec 15 15:30:39 2012

1000000000 digits of pi using Ramanujan formula, 125273397 terms.
        initialization : 15.850 s
      binary splitting : 6177.520 s                
              division : 790.000 s
           square root : 119.790 s
        multiplication : 78.360 s
-----------------------------------
         calculation : 7182.620 s
              output : 2432.270 s

           user time : 9572.514 s
            sys time : 42.387 s
        memory usage : 9204968 KB

Calculation finished at Sat Dec 15 18:10:53 2012

real160m14.581s
user159m32.518s
sys0m42.423s
nitobe@debian64:~/const_new/const-pi-master$ time ./format.sed pi1G-rama.txt >pi1g-rama.txt

real2m41.404s
user2m39.262s
sys0m2.088s
nitobe@debian64:~/const_new/const-pi-master$ time ./const pi-agm -d1000000000 -o pi1G-agm.txt
Calculation started at Sat Dec 15 18:35:39 2012

1000000000 digits of pi using AGM.
        initialization :  0.000 s
           iteration 1 : 178.290 s
           iteration 2 : 278.120 s
           iteration 3 : 278.280 s
           iteration 4 : 278.630 s
           iteration 5 : 277.700 s
           iteration 6 : 279.700 s
           iteration 7 : 279.920 s
           iteration 8 : 277.230 s
           iteration 9 : 280.430 s
          iteration 10 : 281.650 s
          iteration 11 : 280.880 s
          iteration 12 : 280.740 s
          iteration 13 : 280.720 s
          iteration 14 : 283.320 s
          iteration 15 : 280.400 s
          iteration 16 : 280.400 s
          iteration 17 : 281.510 s
          iteration 18 : 280.970 s
          iteration 19 : 280.330 s
          iteration 20 : 280.230 s
          iteration 21 : 280.000 s
          iteration 22 : 280.570 s
          iteration 23 : 280.520 s
          iteration 24 : 278.800 s
          iteration 25 : 276.900 s
          iteration 26 : 275.450 s
          iteration 27 : 278.470 s
          iteration 28 : 291.890 s
          iteration 29 : 255.460 s
          iteration 30 : 225.770 s
              division : 731.840 s
-----------------------------------
         calculation : 8955.190 s
              output : 2350.680 s

           user time : 11110.358 s
            sys time : 195.520 s
        memory usage : 6256556 KB

Calculation finished at Sat Dec 15 21:44:05 2012

real188m25.535s
user185m10.362s
sys3m15.628s
nitobe@debian64:~/const_new/const-pi-master$ time ./format.sed pi1G-agm.txt >pi1g-agm.txt

real2m36.758s
user2m34.898s
sys0m1.868s
nitobe@debian64:~/const_new/const-pi-master$
ご覧のように、Chudnovskyの式、Ramanujanの式、AGMの三種類のアルゴリズムで10億桁を計算してみた。
さて、いよいよ比較なのだが、巨大な一行野郎を、sedで100桁/行フォーマットにした理由が、ここでわかっていただけると思う。diff が使いたかったんだな
nitobe@debian64:~/const_new/const-pi-master$ time diff pi1g.txt pi1g-rama.txt
1,7c1,8
< Calculation started at Sat Dec 15 12:47:03 2012
< 1000000000 digits of pi using Chudnovsky formula, 70513670 terms.
<         initialization : 17.720 s
<       binary splitting : 4126.380 s
<               division : 771.880 s
<            square root : 115.470 s
<         multiplication : 76.110 s
---
> Calculation started at Sat Dec 15 15:33.
> 9 2012
> 1000000000 digits of pi using Ramanujan formula, 125273397 terms.
>         initialization : 15.850 s
>       binary splitting : 6177.520 s
>               division : 790.000 s
>            square root : 119.790 s
>         multiplication : 78.360 s
9c10
<          calculation : 5108.710 s
---
>          calculation : 7182.620 s
10000013,10000017c10000014,10000018
<               output : 2377.740 s
<            user time : 7453.158 s
<             sys time : 33.330 s
<         memory usage : 7567424 KB
< Calculation finished at Sat Dec 15 14:51:49 2012
---
>               output : 2432.270 s
>            user time : 9572.514 s
>             sys time : 42.387 s
>         memory usage : 9204968 KB
> Calculation finished at Sat Dec 15 18:10:53 2012

real0m34.410s
user0m9.657s
sys0m4.880s
nitobe@debian64:~/const_new/const-pi-master$ time diff pi1g.txt pi1g-agm.txt
1,7c1,35
< Calculation started at Sat Dec 15 12:47:03 2012
< 1000000000 digits of pi using Chudnovsky formula, 70513670 terms.
<         initialization : 17.720 s
<       binary splitting : 4126.380 s
<               division : 771.880 s
<            square root : 115.470 s
<         multiplication : 76.110 s
---
> Calculation started at Sat Dec 15 18:35:39 2012
> 1000000000 digits of pi using AGM.
>         initialization :  0.000 s
>            iteration 1 : 178.290 s
>            iteration 2 : 278.120 s
>            iteration 3 : 278.280 s
>            iteration 4 : 278.630 s
>            iteration 5 : 277.700 s
>            iteration 6 : 279.700 s
>            iteration 7 : 279.920 s
>            iteration 8 : 277.230 s
>            iteration 9 : 280.430 s
>           iteration 10 : 281.650 s
>           iteration 11 : 280.880 s
>           iteration 12 : 280.740 s
>           iteration 13 : 280.720 s
>           iteration 14 : 283.320 s
>           iteration 15 : 280.400 s
>           iteration 16 : 280.400 s
>           iteration 17 : 281.510 s
>           iteration 18 : 280.970 s
>           iteration 19 : 283.
> 30 s
>           iteration 20 : 280.230 s
>           iteration 21 : 280.000 s
>           iteration 22 : 280.570 s
>           iteration 23 : 280.520 s
>           iteration 24 : 278.800 s
>           iteration 25 : 276.900 s
>           iteration 26 : 275.450 s
>           iteration 27 : 278.470 s
>           iteration 28 : 291.890 s
>           iteration 29 : 255.460 s
>           iteration 30 : 225.770 s
>               division : 731.840 s
9c37
<          calculation : 5108.710 s
---
>          calculation : 8955.190 s
10000013,10000017c10000041,10000046
<               output : 2377.740 s
<            user time : 7453.158 s
<             sys time : 33.330 s
<         memory usage : 7567424 KB
< Calculation finished at Sat Dec 15 14:51:49 2012
---
>               output : 2350.680 s
>            user time : 11113.
> 58 s
>             sys time : 195.520 s
>         memory usage : 6256556 KB
> Calculation finished at Sat Dec 15 21:44:05 2012

real0m23.249s
user0m9.133s
sys0m3.004s
nitobe@debian64:~/const_new/const-pi-master$ time diff pi1g.txt ../../pi/pi1g.txt
1,9d0
< Calculation started at Sat Dec 15 12:47:03 2012
< 1000000000 digits of pi using Chudnovsky formula, 70513670 terms.
<         initialization : 17.720 s
<       binary splitting : 4126.380 s
<               division : 771.880 s
<            square root : 115.470 s
<         multiplication : 76.110 s
< -----------------------------------
<          calculation : 5108.710 s
10000011,10000017c10000002
< 644
< 
<               output : 2377.740 s
<            user time : 7453.158 s
<             sys time : 33.330 s
<         memory usage : 7567424 KB
< Calculation finished at Sat Dec 15 14:51:49 2012
---
> 64388312

real0m23.689s
user0m9.213s
sys0m3.012s
nitobe@debian64:~/const_new/const-pi-master$
恐るべしdiff。1GByteのファイルを30秒前後で比較している。

参考までにHanhong Xue 氏のファイルと、すずきひろのぶ氏のファイルの頭とケツをお見せしよう。
Chudnovsky formula by Hanhong Xue
       1 Calculation started at Sat Dec 15 12:47:03 2012
       2 1000000000 digits of pi using Chudnovsky formula, 70513670 terms.
       3         initialization : 17.720 s
       4       binary splitting : 4126.380 s
       5               division : 771.880 s
       6            square root : 115.470 s
       7         multiplication : 76.110 s
       8 -----------------------------------
       9          calculation : 5108.710 s
      10 3.
      11 1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679
      12 8214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196
      13 4428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273
     :
     :
10000009 4219776753871319682188195635848934815504410194647387557034502943416861599324354199731814355060392734
10000010 6434543524276655356743570219396394581990548327874671398682093196353628204612755715171395115275045519
10000011 644
10000012 
10000013               output : 2377.740 s
10000014            user time : 7453.158 s
10000015             sys time : 33.330 s
10000016         memory usage : 7567424 KB
10000017 Calculation finished at Sat Dec 15 14:51:49 2012

AGM by Hironobu Suzuki
       1 3.
       2 1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679
       3 8214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196
       4 4428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273
     :
     :
10000000 4219776753871319682188195635848934815504410194647387557034502943416861599324354199731814355060392734
10000001 6434543524276655356743570219396394581990548327874671398682093196353628204612755715171395115275045519
10000002 64388312
xue 氏のファイルは、頭とケツに諸々の情報が入っているので9行ずれている。

4種類の実装で10億2桁目まで同一であることを確認しました。

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

ぎゃぼ!sed スクリプトが、時間の文字列”0.3”で悪さをしている!正しい挙動だ。ま、数値の比較には影響しないので見なかったことにしておこう。

ご希望の方には円周率10億桁のテキストファイルをメールで送付させていただきます。1.01GByteですけど...何か...問題でも?。

 

— posted by nitobe at 06:42 pm   commentComment [0] 

円周率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 Link on BOINCLink 及び、World Community GridLink on BOINC がリソースの75%を占有している。

pi1g.txt
       1 3.
       2 1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679
       3 8214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196

          :
          :
10000000 4219776753871319682188195635848934815504410194647387557034502943416861599324354199731814355060392734
10000001 6434543524276655356743570219396394581990548327874671398682093196353628204612755715171395115275045519
10000002 64388312
さて、この値が正しいかどうか?それが問題だ。つづく

 

— posted by nitobe at 09:35 pm   commentComment [0] 

円周率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上で円周率の計算をおこなうLink 」すずきひろのぶ氏の百万桁計算のソースを使用する。変更点は以下の通り。itchynyさんの修正も適用している。すずきひろのぶ氏が指摘している通り、gmp の precision のサイズの上限により、円周率計算は6.4億桁が上限となる。12/15追記:大嘘です円周率10億桁を計算する #1Link
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は演算繰り返し回数だが、求め得る桁数と、以下の関係がある。
loopmaxcorrect 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

8億桁を指定しているが、6億4千万桁しかない。フォーマット変換後、改行「¥n」が6百万個増えている。百桁/行フォーマットなのでエディタで開ける。重いけど。これは vi で見てみた。
      1 3.
      2 1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679
      3 8214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196
    :
    :
6464569 9336326991629433891751838174109742989418332960270103331518724202832616454999878479262174537890538106
6464570 6455256718445892438033533034196466727204878408211469825653223223178441813002292106453186861390731397
6464571 2309985607826687727041129258325166030790611057861923604585439381176655074474471470989090334958

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

 

— posted by nitobe at 05:34 pm   commentComment [0] 

T: Y: ALL: Online:
ThemeSwitch
  • Basic
Created in 1.1424 sec.
prev
2012.12
next
            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          
 
strawberry-linux geigercounter Ver.2
Sibasaki, Cyofu City, Tokyo, JAPAN
blogBar