digital西行庵blog

pdp11のエミュレータを探してググっていたら、とんでもないものをみつけた。
「Pi with dc 」Julius Schmidt：German natural science student interested in science。unix に標準的に付属しているスタック型（逆ポーランド記法）電卓 dc で実装した spigot algorithm の円周率計算プログラムである。一桁づつメモリのある限り動作し続ける。遅いけど。

echo '1sq180sr60st2si[3li*1+d1+*3*suli27*12-lq*5lr*+lt5*/d48+Psy10lqlid2*1-***10lulqli5*2-*lr+lylt*-**srsqlult*stli1+silmx]smlmx' | dc

このページからもリンクされているが、原典は「Unbounded Spigot Algorithms for the Digits of Pi 」JEREMY GIBBONS：University Lecturer in Software Engineering and Continuing Education at the University of Oxford

> piG3 = g(1,180,60,2) where 
>   g(q,r,t,i) = let (u,y)=(3*(3*i+1)*(3*i+2),div(q*(27*i-12)+5*r)(5*t)) 
>                in y : g(10*q*i*(2*i-1),10*u*(q*(5*i-2)+r-y*t),t*u,i+1)

げげ、Haskell じゃぁないか。

root@jesse:~# apt-get install haskell-platform
パッケージリストを読み込んでいます... 完了
依存関係ツリーを作成しています                
状態情報を読み取っています... 完了
提案パッケージ:
  haskell-platform-doc haskell-platform-prof
以下のパッケージが新たにインストールされます:
  haskell-platform
アップグレード: 0 個、新規インストール: 1 個、削除: 0 個、保留: 3 個。
7,340 B 中 0 B のアーカイブを取得する必要があります。
この操作後に追加で 7,168 B のディスク容量が消費されます。
以前に未選択のパッケージ haskell-platform を選択しています。
(データベースを読み込んでいます ... 現在 187646 個のファイルとディレクトリがインストールされています。)
.../haskell-platform_2013.2.0.0.debian12_all.deb を展開する準備をしています ...
haskell-platform (2013.2.0.0.debian12) を展開しています...
haskell-platform (2013.2.0.0.debian12) を設定しています ...
root@jesse:~# ^dログアウト
hiroaki@jesse:~$ cd pi
hiroaki@jesse:~/pi$ cat piG3.hs
>piG3 = g(1,180,60,2) where
>   g(q,r,t,i) = let (u,y)=(3*(3*i+1)*(3*i+2),div(q*(27*i-12)+5*r)(5*t))
>                in y : g(10*q*i*(2*i-1),10*u*(q*(5*i-2)+r-y*t),t*u,i+1)
main = print piG3
hiroaki@jesse:~/pi$ ghc --make piG3
[1 of 1] Compiling Main             ( piG3.hs, piG3.o )
Linking piG3 ...
hiroaki@jesse:~/pi$ ./piG3
[3,1,4,1,5,9,2,6,5,3,5,8,9,7,9,3,2,3,8,4,6,2,6,4,3,3,8,3,2,7,9,5,0,2,8,8,4,1,9,7,1,6,9,3,9,9,3,7,5,1,0,5,8,2,0,9,7,4,9,4,4,5,9,2,3,0,7,8,1,6,4,0,6,2,8,6,2,0,8,9,9,8,6,2,8,0,3,4,8,2,5,3,^c

げげ、リスト表示だ。こうなったらC言語で実装してやる。

hiroaki@jesse:~/pi$ cat piG3.c
// Pi by spigot algorithm.
//
// [1] JEREMY GIBBONS, 
// "Unbounded Spigot Algorithms for the Digits of Pi"(2005/07/29)
// http://www.cs.ox.ac.uk/people/jeremy.gibbons/publications/spigot.pdf 
// JEREMY GIBBONS is a University Lecturer in Software Engineering and Continuing Education at the University of Oxford.
// > piG3 = g(1,180,60,2) where 
// > g(q,r,t,i) = let (u,y)=(3*(3*i+1)*(3*i+2),div(q*(27*i-12)+5*r)(5*t)) 
// > in y : g(10*q*i*(2*i-1),10*u*(q*(5*i-2)+r-y*t),t*u,i+1)
//
// [2] Julius Schmidt, "AIJU'S HOMEPAGE/Code/Pi"
// http://aiju.de/code/pi 
// Pi with dc
// echo '1sq180sr60st2si[3li*1+d1+*3*suli27*12-lq*5lr*+lt5*/d48+Psy10lqlid2*1-***10lulqli5*2-*lr+lylt*-**srsqlult*stli1+silmx]smlmx' | dc
// It's a great implementation.
// 
// Implementation of C with GMP by Hiroaki Nitobe.　2016/02/07
//

#include        <stdio.h>
#include        <signal.h>
#include        <gmp.h>
#include        "time_tools.h"

static int detect_sigint = 0;

void signal_handler(int signo){
        if(signo==SIGINT){
        fprintf(stderr, "\nReceived SIGINT\n");
        detect_sigint = 1;
}
}

void g(void){
        mpz_t        q, r, t, y, temp;
        unsigned long int        i, u, d;

        mpz_init(q);
        mpz_init(r);
        mpz_init(t);
        mpz_init(y);
        mpz_init(temp);

        // q=1
        mpz_set_ui(q, 1UL);
        // r=180
        mpz_set_ui(r, 180UL);
        // t=60
        mpz_set_ui(t, 60UL);
        // i=2
        i = 2UL;
        // d=0
        d = 0UL;

        while(1) {
        // u=3*(3*i+1)*(3*i+2)
                u=3UL*(3UL*i+1UL)*(3UL*i+2UL);
        // y=(q*(27*i-12)+5*r)/(5*t)
                mpz_mul_ui(y, q, 27UL*i-12UL);
                mpz_mul_ui(temp, r, 5UL); 
                mpz_add(y, y, temp);
                mpz_mul_ui(temp, t, 5UL);
                mpz_fdiv_q(y, y, temp);
                mpz_out_str(stdout, 10, y);
                                if(d%100UL==0UL){
                                        if(d==0){
                                                printf(".\n");
                                        }else{
                                                printf("\n");
                                        }
                                        if(d%10000UL==0UL){
                                                fprintf(stderr, "%lu, ", d);
                                                my_timestamp(stderr);
                                                if(detect_sigint==1)break;
                                        }
                                }

        // next r = 10*u*(q*(5*i-2)+r-y*t)
                mpz_mul_ui(temp, q, 5UL*i-2UL);
                mpz_add(r, r, temp);
                mpz_mul(temp, y, t);
                mpz_sub(r, r, temp);
                mpz_mul_ui(r, r, 10UL*u);
        // next q = 10*q*i*(2*i-1)
                mpz_mul_ui(q, q, 10UL*i*(2UL*i-1));
        // next t = t*u
                mpz_mul_ui(t, t, u);
        // next i = i+1
                i = i + 1;
                d = d + 1;
        }
        mpz_clear(q);
        mpz_clear(r);
        mpz_clear(t);
        mpz_clear(y);
        mpz_clear(temp);
        fprintf(stderr, "Normal end.\n");
}

int main(int argc, char *argv[]){
        if(signal(SIGINT, signal_handler) == SIG_ERR){
                fprintf(stderr, "Can't catch SIGINT\n");
        }

        fprintf(stderr, "with gcc-%d.%d.%d, gmp-%s\n",
                        __GNUC__, __GNUC_MINOR__, __GNUC_PATCHLEVEL__, gmp_version);
                                my_date(stderr);
                                my_set_start(stderr);

        g();
                                my_date(stderr);
        return 0;
}
hiroaki@jesse:~/pi$ 
hiroaki@jesse:~/pi$ cat time_tools.c
#include        <stdio.h>
#include        <stdlib.h>
#include        <time.h>
#include        <sys/time.h>
#include        <sys/types.h>
#include        <unistd.h>

void my_date(FILE *fp);
double my_get_time(void);
void my_set_start(FILE *fp);
void my_timestamp(FILE *fp);

static double        start, lap;

void my_date(FILE *fp) {
        time_t        timer;
        time(&timer);
        fprintf(fp, "%s", ctime(&timer));
}

double my_get_time(void) {
        struct timeval        tv;

        gettimeofday(&tv, NULL);

        return (tv.tv_sec + tv.tv_usec * 1e-6);
}

void my_set_start(FILE *fp) {
        start = lap = my_get_time();
        fprintf(fp, "split: %f, lap: %f\n", start - start, lap - lap);
}

void my_timestamp(FILE *fp) {
        double        now;

        now = my_get_time();

        fprintf(fp, "split: %f, lap: %f\n", now - start, now - lap);

        lap = now;
}
hiroaki@jesse:~/pi$ 
hiroaki@jesse:~/pi$ cat time_tools.h
extern void my_date(FILE *fp);
extern double my_get_time(void);
extern void my_set_start(FILE *fp);
extern void my_timestamp(FILE *fp);
hiroaki@jesse:~/pi$ 
hiroaki@jesse:~/pi$ gcc -l gmp -o piG3 piG3.c time_tools.c
hiroaki@jesse:~/pi$ ./piG3 > piG.txt
with gcc-4.9.2, gmp-6.0.0
Thu Feb 11 17:50:18 2016
split: 0.000000, lap: 0.000000
0, split: 0.000048, lap: 0.000048
10000, split: 0.260359, lap: 0.260311
20000, split: 1.232651, lap: 0.972292
30000, split: 3.073861, lap: 1.841210
40000, split: 6.088222, lap: 3.014361
^C
Received SIGINT
50000, split: 10.240260, lap: 4.152038
Normal end.
Thu Feb 11 17:50:28 2016
hiroaki@jesse:~/pi$ 

んん。動いてるみたいだ。
さて、最初の dc 版をみてみよう。
　
echo '1sq180sr60st2si[3li*1+d1+*3*suli27*12-lq*5lr*+lt5*/d48+Psy10lqlid2*1-***10lulqli5*2-*lr+lylt*-**srsqlult*stli1+silmx]smlmx' | dc

1sq	１をｑレジスタに積む	g(1,180,60,2)
180sr	１８０をｒレジスタに積む	g(1,180,60,2)
60st	６０をｔレジスタに積む	g(1,180,60,2)
2si	２をｉレジスタに積む	g(1,180,60,2)
[	ここから
3li*	3 とｉレジスタを掛けてメインスタックに積む	3*(3*i+1)*(3*i+2)
1+	メインスタックに１を足す	3*(3*i+1)*(3*i+2)
d	メインスタックを複製して積む	3*(3*i+1)*(3*i+1+1)
1+	メインスタックに１足す	3*(3*i+1)*(3*i+1+1)
*	メインスタック（top と next）を掛ける	3*(3*i+1)*(3*i+2)
3*	メインスタックに３を掛ける	3*(3*i+1)*(3*i+2)
su	メインスタックをｕレジスタに積む	let (u,y)=
li27*	ｉレジスタと２７を掛けてメインスタックに積む	div(q*(27*i-12)+5*r)(5*t)
12-	メインスタックから１２引く	div(q*(27*i-12)+5*r)(5*t)
lq*	メインスタックにｑレジスタを掛ける	div(q*(27*i-12)+5*r)(5*t)
5lr*	５とｒレジスタを掛けメインスタックに積む	div(q*(27*i-12)+5*r)(5*t)
+		div(q*(27*i-12)+5*r)(5*t)
lt5*		div(q*(27*i-12)+5*r)(5*t)
/		div(q*(27*i-12)+5*r)(5*t)
d	印刷用にメインスタックを複製
48+	３０ｈ（４８ｄ）の下駄を履かせて	asciiコード変換
P	プリント	y :
sy	メインスタックをｙレジスタに積む	let (u,y)=
10lqlid2*1-***	１０とｑとｉをメインスタックに積んでｉを複製して２を掛けて１を引いて掛けて掛けて掛ける	g(10*q*i*(2*i-1),10*u*(q*(5*i-2)+r-y*t),t*u,i+1)
10lulqli5*2-*lr+lylt*-**	１０とｕとｑとｉと５をメインスタックに積んで２を引いて掛けてｒを足してｙとｔを掛けて引いて掛けて掛ける	g(10*q*i*(2*i-1),10*u*(q*(5*i-2)+r-y*t),t*u,i+1)
srsq	メインスタックのtopをｒレジスタに、nextをｑレジスタに積む	g(10*q*i*(2*i-1),10*u*(q*(5*i-2)+r-y*t),t*u,i+1)
lult*st	ｕレジスタとｔレジスタを掛けてｔレジスタに積む	t*u
li1+si	ｉレジスタと１を足してｉレジスタに積む	i+1
lmx	ｍレジスタのマクロを無条件で実行	再帰マクロ呼び出し
]sm	ここまでをマクロとしてｍレジスタに積む	マクロ定義
lmx	ｍレジスタのマクロを無条件で実行	初回マクロ呼び出し

以上の文字列をechoしてdcに食わせるという寸法だ。
It's a great implementation.

pi  

実は debian には pi コマンドがある。

root@Wheezy64:~# apt-get install pi
パッケージリストを読み込んでいます... 完了
　　　　：
　　【中略】
　　　　：
pi (1.3.2-1.2) を設定しています ...
root@Wheezy64:~# time pi 1000000000 > pi_1g.txt

real94m3.688s
user92m36.915s
sys1m20.685s
root@Wheezy64:~#

インストールして桁数指定して実行するだけ。10億桁を5643.688secで計算しやがった。(注：例によって裏では宇宙人探しと癌DNA解析とエイズ分子構造計算をやっているので計算時間はあてにならない。）gmpではなくclnライブラリを使っているらしい。ソースは、/usr/src/cln-1.3.2/src/float/transcendentalにあるので、興味がある人は見てみるといい。cl_LF_pi.cc。

pi  

---

まったく意味がないことを承知で捏ね繰り廻す。何をやっているかは、シェルスクリプトを解読し給へ。

nitobe@debian64:~/itchyny$ cat count1.sh
#!/bin/sh
start=1000001
step=1000000
end=10000001
for i in `seq $start $step $end`
do
        echo -n "$i, "
        for k in `seq 0 9`
        do
                sed -n -e "2, $i p" pi1g.txt| tr -dc $k| wc -c| tr '¥n' ','
        done
        echo
done
nitobe@debian64:~/itchyny$ chmod +x count1.sh
nitobe@debian64:~/itchyny$ time ./count1.sh
1000001, 9999922,10002475,10001092,9998442,10003863,9993478,9999417,9999610,10002180,9999521,
2000001, 19997437,20003774,20002185,20001410,19999846,19993031,19999161,20000287,20002307,20000562,
3000001, 29998356,30000582,30006337,29999867,29999810,29993099,29998913,29999071,30003683,30000282,
4000001, 39996048,39997375,40011791,39995030,40001014,39992123,40001899,40000314,40005735,39998671,
5000001, 49995279,50000437,50011436,49992409,50005121,49990678,49998820,50000320,50006632,49998868,
6000001, 59991725,59997597,60008591,59992558,60007991,59990211,60003895,59998772,60010958,59997702,
7000001, 69989891,69997755,70006497,69994028,70009581,69994537,70003795,69997014,70005161,70001741,
8000001, 79991897,79997003,80003316,79989651,80016073,79996120,80004148,79995109,80002933,80003750,
9000001, 89991208,89998381,90000968,89990083,90013132,89996086,90006412,89995658,90001979,90006093,
10000001, 99993942,99997334,100002410,99986911,100011958,99998885,100010387,99996061,100001839,100000273,

real4m12.549s
user5m20.528s
sys1m45.579s
nitobe@debian64:~/itchyny$

一億桁づつの数字の累積出現回数だぜぃ。累積度数分布というのかな？

nitobe@debian64:~/itchyny$ cat count2.sh
#!/bin/sh
size=1000000
for i in `seq 2 $size 10000001`
do
        j=$((i+size-1))
        echo -n "$i, $j, "
        for k in `seq 0 9`
        do
                sed -n -e "$i,$j p" pi1g.txt| tr -dc $k| wc -c| tr '¥n' ','
        done
        echo
done
nitobe@debian64:~/itchyny$ chmod +x count2.sh
nitobe@debian64:~/itchyny$ time ./count2.sh
2, 1000001, 9999922,10002475,10001092,9998442,10003863,9993478,9999417,9999610,10002180,9999521,
1000002, 2000001, 9997515,10001299,10001093,10002968,9995983,9999553,9999744,10000677,10000127,10001041,
2000002, 3000001, 10000919,9996808,10004152,9998457,9999964,10000068,9999752,9998784,10001376,9999720,
3000002, 4000001, 9997692,9996793,10005454,9995163,10001204,9999024,10002986,10001243,10002052,9998389,
4000002, 5000001, 9999231,10003062,9999645,9997379,10004107,9998555,9996921,10000006,10000897,10000197,
5000002, 6000001, 9996446,9997160,9997155,10000149,10002870,9999533,10005075,9998452,10004326,9998834,
6000002, 7000001, 9998166,10000158,9997906,10001470,10001590,10004326,9999900,9998242,9994203,10004039,
7000002, 8000001, 10002006,9999248,9996819,9995623,10006492,10001583,10000353,9998095,9997772,10002009,
8000002, 9000001, 9999311,10001378,9997652,10000432,9997059,9999966,10002264,10000549,9999046,10002343,
9000002, 10000001, 10002734,9998953,10001442,9996828,9998826,10002799,10003975,10000403,9999860,9994180,

real2m49.034s
user2m44.614s
sys0m36.742s
nitobe@debian64:~/itchyny$

こっちは一億桁づつの数字の区間出現回数だぜぃ。区間度数分布？

ちっともワイルドじゃない。

コピペして、Windows上のテキストファイルにして、excelで取り込んでグラフ化してみた。もうひと捻りほしいところだ。宿題。

---

2013/01/10追加：平均値との差

---

---

---

注：３．１４・・・の最初の３はカウントしていないから、１足してね。

---

xlsx ファイル置いときます。好きにして頂戴。2013/01/10 count1.xlsx差し替え。

pi  

---

---

2012-12-18 付で、多倍長演算ライブラリ gmp が、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.bz2 
     :
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# 

まずは、円周率１億桁の実行時間を比較する。使用するソースは itchyny さんのもの。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万桁のベンチマークは、「円周率ベンチマークをやり直す 」を参照していただきたい。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 不能。

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

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

pi  

---

---

と、いうわけで、検証。
macroxue / const-pi のzipファイルを入手する。
件の Hanhong Xue 氏の円周率計算プログラム、新しめのバージョンである。

nitobe@debian64:~$ wget https://github.com/macroxue/const-pi/archive/master.zip 
     :
     :
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億桁の計算はできるのだが、２進１０進変換で落ちる。たぶんメモリ不足だ。フォーマットも非常に凝ったものだが、扱いにくいので、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の三種類のアルゴリズムで１０億桁を計算してみた。
さて、いよいよ比較なのだが、巨大な一行野郎を、sedで１００桁／行フォーマットにした理由が、ここでわかっていただけると思う。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のファイルを３０秒前後で比較している。

参考までに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 氏のファイルは、頭とケツに諸々の情報が入っているので９行ずれている。

４種類の実装で１０億２桁目まで同一であることを確認しました。

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

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

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

pi