■BBPとは？
BBPは円周率の小数点以下の任意の桁の値を求める事が出来るアルゴリズムです。
参考：BBPアルゴリズム
http://www14.ocn.ne.jp/~kk62526/pi/BBP.html
参考：円周率の望みの桁を直接計算する法
http://d.hatena.ne.jp/rikunora/20110627/p1
参考：円周率の任意の桁を一発で計算する
http://www.riverplus.net/sci/2005/05/post_83.html
参考：π計算プログラム 「スーパーπ Ver 1.1」
http://www1.coralnet.or.jp/kusuto/PI/super_pi.html
■以下にCLISPのコードがあった。
※BBPのメリットを無視するようだが、
単純に円周率を毎回1から計算する方針で進める。
参考：GNU/Linux上で円周率の計算をおこなう
http://h2np.net/pi/
■しかしCLISPには組み込みの「pi」関数がある
$ echo "pi" | clisp -q | tail -1
3.1415926535897932385L0
■というわけでmypiに変更
(setf (EXT:LONG-FLOAT-DIGITS) 330)
(defun mypi ()
(defun s(k)
(* (/ 1.0L0 (expt 16.0L0 k))
(- (/ 4.0L0 (+ 1.0L0(* 8.0L0 k)))
(+ (/ 2.0L0 (+ 4.0L0 (* 8.0L0 k)))
(+ (/ 1.0L0 (+ 5.0L0 (* 8.0L0 k)))
(/ 1.0L0 (+ 6.0L0 (* 8.0L0 k))))))))
(defun m(n)
(if (> n 90L0) 0.0L0
(+ (m (+ n 1.0L0))
(s n))))
(m 0.0L0))
■時間を計る時には、「-o」オプションが使えなかったので、リダイレクトする。
$ time (echo "(mypi)" | clisp -q -i bbp.lisp > bbp.log)
real 0m0.115s
user 0m0.008s
sys 0m0.040s
■小数点以下100桁を取り出します。
$ tail -1 bbp.log | awk -F\. '{print $2}' | awk -F L '{print $1}' | wc -c
107
$ echo "(mypi)" | clisp -q -i bbp.lisp | tail -1 | \
cut -c 3-102 | sed s/"."/"&\n"/g | \
awk '{if (NR % 10 == 0) print $1 ":"; else print $1}' | \
xargs echo -n | sed s/" "//g | sed s/":"/"\n"/g;echo ""
1415926535
8979323846
2643383279
5028841971
6939937510
5820974944
5923078164
0628620899
8628034825
3421170679
■有効桁数については、「EXT:LONG-FLOAT-DIGITS」の
フローティングポイントを変更すればよいので、以下のようにします。
$ cat bbp.lisp | sed s/"330"/"1000"/ > bbp1.lisp
$ cat bbp.lisp | sed s/"330"/"10000"/ > bbp2.lisp
$ cat bbp.lisp | sed s/"330"/"100000"/ > bbp3.lisp
$ cat bbp.lisp | sed s/"330"/"1000000"/ > bbp4.lisp
$ cat bbp.lisp | sed s/"330"/"3300000"/ > bbp5.lisp
■毎回以下のようにするのは面倒。
$ cat bbp3.log | awk -F\. '{print $2}' | awk -F\L '{print $1}' | tail -1 | wc -c
309
■以下で充分。
$ tail -1 bbp1.log | wc -c | awk '{print $1-($1%100)}'
300
■ループでbbp3.lispまで計算してみる。
$ for list in bbp[1-3].lisp;do \
LOG=$(echo $list | sed s/"\.lisp"/"\.log"/); \
time (echo "(mypi)" | clisp -q -i $list > $LOG); \
done
real 0m0.155s
user 0m0.012s
sys 0m0.052s
real 0m1.183s
user 0m0.308s
sys 0m0.124s
real 3m7.682s
user 1m5.108s
sys 0m2.260s
■bbp4.lispは約2日かかった。bbp5.lispは現在も実行中。
$ ps -ef | grep bbp5.lisp | grep -v grep | awk '{print $7}'
3-23:41:21
■有効桁数を見てみる。
※最近のCentOSでは、「tail -1」が通らない。
その場合は「tail -n 1」とする必要がある。
$ for list in bbp[1-3].log;do \
tail -1 "$list" | wc -c | awk '{print $1-($1%100)}'; \
done
300
3000
30100
■bbp4.logの有効桁数は以下の通り、30万桁。
$ tail -1 bbp5.log | wc -c | awk '{print $1-($1%100)}'
301000
■bbp5.logは、100万桁の予定。
以下の通りCPU使用率は高いがメモリは余り使っていない。
$ top -b -n 1 | head -7 | tail -1;top -b -n 1 | grep lisp
PID USER PR NI VIRT RES SHR S %CPU %MEM TIME+ COMMAND
9183 labunix 20 0 113m 69m 2424 R 99.9 13.8 5749:54 lisp.run
■以下の通りメモリ512MBしか割り当てていない仮想マシン
$ free
total used free shared buffers cached
Mem: 511004 452172 58832 0 118972 109420
-/+ buffers/cache: 223780 287224
Swap: 1421712 8308 1413404
■以下の通りに変換。
※「５．」で改行とスペースを消すので、
それぞれ別の記号に前変換してエスケープしている。
１．最終行を取り出し
２．「3.1」を除いた30万1千桁の取り出し、すべての文字の後ろに改行を挿入⇒列にする
３．100列ごとに「%」を挿入、それ以外の10列ごとに「:」を挿入
５．「\n」を取り除き、余分なスペースを削除
６．「:」をスペースに、「%」を改行に変換
$ tail -1 bbp4.log | \
cut -c 3-301002 | sed s/"."/"&\n"/g | \
awk '{if (NR % 100 == 0) print $1 "%"; else if (NR % 10 == 0) print $1 ":"; else print $1}' | \
xargs echo -n | sed s/" "//g | \
sed s/":"/" "/g | sed s/"%"/"\n"/g > bbp4.txt
■10桁区切り、10回が一行で、「wc -l」で行数を掛ける
$ wc -l bbp4.txt | awk '{print 10*10*$1}'
301000
■ワード単位(英文のスペース区切りの意味)で1ワード=(10桁)として数える。
$ wc -w bbp4.txt | awk '{print 10*$1}'
301000
■サイズは328k
$ du -h bbp4.txt
328K bbp4.txt
■載せようとしたら文字数オーバーでした。。。
$ cat bbp4.log
■計算したpiをチェックする
※参考1のダウンロードサイトはサービスが終了しておりダウンロード出来ない。
参考1：π1000万桁のDownload
http://www1.coralnet.or.jp/kusuto/PI/PI-data/index.html
参考2：円周率 1,000,000 桁
http://www.kisaragiweb.jp/pi/pi1m.htm
■チェックしたが、違うらしい。。。
$ for num in `seq 1 4`;do \
echo -n "bbp$num :" ; \
diff bbp${num}.txt check${num}.txt | head -1; \
done
bbp1 :4d3
bbp2 :31d30
bbp3 :302d301
bbp4 :2,3010c2,3010
■チェックファイルは以下のようにして保存。
※一行で出来る。
$ w3m -dump -cols 110 http://www.kisaragiweb.jp/pi/pi1m.htm | \
grep "^[0-9]" | \
sed s/" "/"\n"/g > check.txt
■100万桁あることを確認。
$ wc -l check.txt | awk '{print 10*$1}'
1000000
■10桁づつ出力
$ for list in bbp[1-4].log;do
TXT=$(echo $list | sed s/"\.log"/"\.txt"/); \
tail -n 1 $list | \
awk -F\. '{print $2}' | awk -F L '{print $1}' | sed s/"."/"&\n"/g | \
awk '{if (NR % 10 == 0) print $1 ":"; else print $1}' | \
xargs echo -n | sed s/" "//g | \
sed s/":"/"\n"/g > $TXT
done
■check用のファイル切り出し。
$ for list in bbp*.txt;do \
LINE=$(wc -l $list | awk '{print $1}') ; \
CHECKLOG=$(echo $list | sed s/"\.txt"/"\.check"/); \
head -n $LINE check.txt > $CHECKLOG; \
done
■lisp自体の再実行に時間をかけるのもアレなので、コンパイルしよう。
$ for list in bbp[1-4].lisp;do clisp -q -c "$list";done | head -3
;; Compiling file ~/bbp1.lisp ...
;; Wrote file ~/bbp1.fas
0 errors, 0 warnings
■ついでに性能の高いシステム側で再実行
※仮想マシンで2日かかった処理が8分で終わった。。。
$ for list in bbp[1-4].fas;do \
LOG=$(echo $list | sed s/"\.fas"/"\.log"/); \
echo "(mypi)" | clisp -q -i "$list" > $LOG; \
done
■再確認したが、やはり違うようだ。
$ for num in `seq 1 4`;do \
echo -n "bbp${num} : "; \
head -n 1 "bbp${num}.lisp" | awk '{print $3}' | sed s/")"//g | xargs echo -n; \
echo -n " : "; \
diff bbp${num}.txt check${num}.txt | head -1; \
done
bbp1 : 1000 : 4d3
bbp2 : 10000 : 31d30
bbp3 : 100000 : 302d301
bbp4 : 1000000 : 2,3010c2,3010
■100桁まで有効と書いてある。。。
参考：clispで円周率計算のプログラムを書いてみた
http://h2np.net/tips/wiki/index.php?clisp%A4%C7%B1%DF%BC%FE%CE%A8
■なるほど、確かに、すべて同じ間違いをしています。
$ for num in `seq 1 4`;do \
diff -y -W 80 bbp${num}.txt bbp${num}.check | \
nl -n rz | grep "|" | head -n 1; \
done
000012 3281510877 | 3282306647
000012 3281510877 | 3282306647
000012 3281510877 | 3282306647
■まとめ。
今回は時間の無駄だったということで。。。
