■wikiのcommon lispプログラムをclispで試す
http://ja.wikipedia.org/wiki/%E7%B4%A0%E6%95%B0%E5%88%A4%E5%AE%9A
■素数判定(エラトステネスの篩)
プログラム部分をコピペ
$ cat prime.lisp
(defun eratosthenes (n)
(labels ((foo (i ls)
(if (= i 1)
ls
(foo (1- i) (cons i ls))))
(spam (ls0 ls1)
(if (> (expt (car ls1) 2) (car (last ls0)))
(revappend ls1 ls0)
(let ((ls2
(remove-if #'(lambda (x)
(zerop (mod x (car ls1))))
ls0)))
(spam (cdr ls2) (cons (car ls2) ls1))))))
(let ((lst (foo n nil)))
(spam (cdr lst) (list (car lst))))))
■timeコマンドに渡す
$ time echo "(eratosthenes 1000)" | clisp -q -i prime.lisp
;; Loading file prime.lisp ...
;; Loaded file prime.lisp
[1]>
(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103
107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199
211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307 311 313
317 331 337 347 349 353 359 367 373 379 383 389 397 401 409 419 421 431 433
439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541 547 557 563
569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659 661 673
677 683 691 701 709 719 727 733 739 743 751 757 761 769 773 787 797 809 811
821 823 827 829 839 853 857 859 863 877 881 883 887 907 911 919 929 937 941
947 953 967 971 977 983 991 997)
real 0m0.077s
user 0m0.012s
sys 0m0.056s
■以下「()」内は10引いた数を記載する
$ echo ";; Loading file prime.fas ...
> ;; Loaded file prime.fas
> [1]>
> " | wc -w
10
■10までの素数の数
$ time echo "(eratosthenes 10)" | clisp -q -i prime.lisp | wc -w
14(4)
real 0m1.400s
user 0m0.016s
sys 0m0.348s
■100までの素数の数
$ time echo "(eratosthenes 100)" | clisp -q -i prime.lisp | wc -w
35(25)
real 0m1.375s
user 0m0.024s
sys 0m0.232s
■1000までの素数の数
$ time echo "(eratosthenes 1000)" | clisp -q -i prime.lisp | wc -w
178(168)
real 0m0.068s
user 0m0.024s
sys 0m0.044s
■オーバーフローのタイミングを確認する
$ time for num in `seq 6500 1 6700`;do \
echo "(eratosthenes $num)" | clisp -q -i prime.lisp 2>&1 | \
grep overflow && echo "$num" && break; \
done
*** - Program stack overflow. RESET
6533
real 0m3.640s
user 0m1.876s
sys 0m1.592s
⇒6532までは正常に表示できます。
■6532までの素数は854です。
$ time echo "(eratosthenes 6532)" | clisp -q -i prime.lisp | wc -w
854
real 0m0.112s
user 0m0.056s
sys 0m0.052s
■コンパイルと最適化について
$ man clisp | grep "\-c lisp-file\$" -A 3
-c lisp-file
Compiles the specified lisp-files to bytecode (*.fas). The compiled
files can then be LOAD[30]ed instead of the sources to gain
efficiency.
■prime.lispをコンパイルし、prime.fasを生成する
$ clisp -c prime.lisp
■エラーが出なくなった
$ time for num in `seq 6500 1 6700`;do \
echo "(eratosthenes $num)" | clisp -q -i prime.fas 2>&1 | \
grep overflow && echo "$num" && break; \
done
real 0m12.621s
user 0m5.348s
sys 0m7.136s
■以下まで数えられるだろうか
参考：10億までの素数の数
http://metalroyd.blog101.fc2.com/blog-date-201108.html
■10000までの素数の数
$ time echo "(eratosthenes 10000)" | clisp -q -i prime.fas | wc -w
1239(1229)
real 0m0.147s
user 0m0.024s
sys 0m0.056s
■100000までの素数の数
$ time echo "(eratosthenes 100000)" | clisp -q -i prime.fas | wc -w
9602(9592)
real 0m0.646s
user 0m0.452s
sys 0m0.188s
■1百万までの素数の数
$ time echo "(eratosthenes 1000000)" | clisp -q -i prime.fas | wc -w
78508(78498)
real 0m10.450s
user 0m6.988s
sys 0m2.228s
■1千万までの素数の数
急激に時間が必要になります。
$ time echo "(eratosthenes 10000000)" | \
clisp -q -i prime.fas > 1sen-man.out 2> 1sen-man.err
real 3m6.092s
user 2m21.433s
sys 0m22.589s
「-o」で指定してもエラーが出ますが、リダイレクトしてみると、
実行中に別の端末から「tail -f」してみると、きちんと動作しているようです。
$ tail -f 1sen-man.out
;; Loading file prime.fas ...
;; Loaded file prime.fas
[1]>
$ wc -w 1sen-man.out
664589(664579)
⇒実行結果に差が出ました。
■再確認
参考：素数定理
http://ja.wikipedia.org/wiki/%E7%B4%A0%E6%95%B0%E5%AE%9A%E7%90%86
⇒こちらと一致しました。
■1億までの素数の数
だいぶかかりそうです。
$ tail -f 1oku.out
;; Loading file prime.fas ...
;; Loaded file prime.fas
[1]> ^C
■作業中の経過時間の確認
$ ps -ef | grep clisp | cut -c 40- | awk '{print $1 " " $2}'
00:07:15 /usr/lib/clisp-2.44.1/base/lisp.run
00:00:00 grep
予想だと90分ほどかかると思います。
さて、5761455に一致するかどうか。
結果が出たら追記します。
--------------------------------------------------------------------------------
□追記
1億は強制終了となりました。。。
$ time echo "(eratosthenes 100000000)" | \
clisp -q -i prime.fas > 1oku.out 2> 1oku.err
強制終了
real 62m2.896s
user 0m57.516s
sys 10m34.324s
□やり方を変える1
範囲の素数リストを小分けにする
参考：(題名無し)
http://ideone.com/LEG4V
□やり方を変える2
参考：10兆までの素数リスト
http://d.hatena.ne.jp/machy3/20110111
□強制終了の理由は負荷だったようだ。
※別の性能が高いマシンに変更
$ time echo "(eratosthenes 100000000)" | \
clisp -q -i prime.fas > prime.out 2> prime.err
real 34m5.555s
user 31m8.877s
sys 1m34.378s
$ wc -w prime.out | awk '{print $1-10}'
5761455
□メモリには余裕があるが、CPU使用率が100%のままだった。
・性能が低いマシン(amd64x2だが、仮想マシン内なので、1コア認識)
$ cat /proc/cpuinfo | grep "MHz\|bogomips\|processor"
processor : 0
cpu MHz : 2411.285
bogomips : 4828.05
・性能が高いマシン(PentiumDualCoreホストsqueeze、2コア認識)
$ cat /proc/cpuinfo | grep "MHz\|bogomips\|processor"
processor : 0
cpu MHz : 1203.000
bogomips : 5185.87
processor : 1
cpu MHz : 1203.000
bogomips : 5186.94
⇒どちらにしろ、少し別のアプローチも考えなければ。。。
□変わったアプローチ
CERN Common Lisp Users Group
http://lisp.cern.ch/guide.html
Common Lisp vs. Scheme macros
http://eli.thegreenplace.net/2007/09/16/common-lisp-vs-scheme-macros/
⇒CERNのアプローチが面白そうだ。
□CERONのコードは、clispの関数単位のコンパイルが盛り込まれている。
$ cat ceron_prime.sh
#!/usr/bin/clisp
(defun primep (x)
"Predicate to test the primality of x."
(let ((sqrt-x (sqrt x)))
(do ((i 2 (1+ i)))
((> i sqrt-x) t)
(when (eq (mod x i) 0)
(return nil)))))
(defun find-primes (m n)
"Find and print n prime numbers greater than m."
(do ((i (1+ m) (1+ i))
(found 0))
((>= found n) t)
(when (primep i)
(format t "~A~%" i)
(incf found))))
(defun print-usage-info ()
"Print usage hints for this script."
(format t "Usage: primes <M> <N>~%")
(format t "Description: Print first N prime numbers greater than M.~%")
(format t "Note: Both arguments M and N must be integers.~%"))
(defun main (args)
"Main function that analyzes command-line arguments
and takes appropriate action."
(if (null (second args))
(print-usage-info)
(let ((arg1 (parse-integer (first args) :junk-allowed t))
(arg2 (parse-integer (second args) :junk-allowed t)))
(if (and (integerp arg1)
(integerp arg2))
(find-primes arg1 arg2)
(print-usage-info)))))
(compile 'find-primes)
(compile 'primep)
(main *args*)
□1から素数4つのリスト
$ time ./ceron_prime.sh 1 4 | column > 1.out
real 0m0.020s
user 0m0.008s
sys 0m0.012s
□10から素数(25-4)つのリスト
$ time ./ceron_prime.sh 10 21 | column > 10.out
real 0m0.029s
user 0m0.024s
sys 0m0.008s
□100から素数(168-25)つのリスト
$ time ./ceron_prime.sh 100 143 | column > 100.out
real 0m0.022s
user 0m0.016s
sys 0m0.008s
□1000から(1229-168)つのリスト
$ time ./ceron_prime.sh 1000 `echo "1229" | awk '{print $1-168}'` | \
column > 1000.out
real 0m0.052s
user 0m0.032s
sys 0m0.024s
□10000から(9592-1229)つのリスト
$ time ./ceron_prime.sh 10000 `echo "9592" | awk '{print $1-1229}'` | \
column > 10000.out
real 0m0.486s
user 0m0.440s
sys 0m0.052s
□100000から(78498-9592)つのリスト
$ time ./ceron_prime.sh 100000 `echo "78498" | awk '{print $1-9592}'` | \
column > 100000.out
real 0m9.245s
user 0m8.749s
sys 0m0.476s
□1000000から(664579-78498)つのリスト
$ time ./ceron_prime.sh 1000000 `echo "664579" | awk '{print $1-78498}'` | \
column > 1000000.out
real 3m24.091s
user 3m23.061s
sys 0m2.556s
□10000000から(5761455-664579)つのリスト
※ここまではコンパイル済みの「prime.fas」とほぼ同時間の為、
最大約32分かかるはずと予想したが、
既に40分が経過した。
□バイト数は、ceron_prime.fasの方が大きい。
$ ls -l prime.fas ceron_prime.fas | cut -c 29-
4254 2011-12-04 03:55 ceron_prime.fas
1723 2011-12-04 03:41 prime.fas
□出力ファイルのバイト数を見るのも面白い
※1億までの処理は実行中
$ wc -w 1*.out
4 1.out
21 10.out
143 100.out
1061 1000.out
8363 10000.out
68906 100000.out
586081 1000000.out
0 10000000.out
664579 合計
□各桁までの素数の数は、上から足したものだ。
各桁なので、足す数=「0の数」となる。
$ wc -w 1*.out | head -2 | awk '{print $1}' | xargs echo -n | \
awk '{print $1+$2}'
25
$ wc -w 1*.out | head -3 | awk '{print $1}' | xargs echo -n | \
awk '{print $1+$2+$3}'
168
$ wc -w 1*.out | head -4 | awk '{print $1}' | xargs echo -n | \
awk '{print $1+$2+$3+$4}'
1229
$ wc -w 1*.out | head -5 | awk '{print $1}' | xargs echo -n | \
awk '{print $1+$2+$3+$4+$5}'
9592
$ wc -w 1*.out | head -6 | awk '{print $1}' | xargs echo -n | \
awk '{print $1+$2+$3+$4+$5+$6}'
78498
□前の桁から見た素数の数の比率を見ると、約0.5づつ増加しているように見える。
$ wc -w 1*.out | head -2 | awk '{print $1}' | xargs echo -n | \
awk '{print ($1+$2)/($1)}'
6.25
$ wc -w 1*.out | head -3 | awk '{print $1}' | xargs echo -n | \
awk '{print ($1+$2+$3)/($1+$2)}'
6.72
$ wc -w 1*.out | head -4 | awk '{print $1}' | xargs echo -n | \
awk '{print ($1+$2+$3+$4)/($1+$2+$3)}'
7.31548
$ wc -w 1*.out | head -5 | awk '{print $1}' | xargs echo -n | \
awk '{print ($1+$2+$3+$4+$5)/($1+$2+$3+$4)}'
7.80472
$ wc -w 1*.out | head -6 | awk '{print $1}' | xargs echo -n | \
awk '{print ($1+$2+$3+$4+$5+$6)/($1+$2+$3+$4+$5)}'
8.18369
□次の比率を8.5とした場合
$ echo "78498" | awk '{print $1*8.5}'
667233
□以下の表を見ると、7桁で664579つの素数がある。
Prime Counting Function
http://mathworld.wolfram.com/PrimeCountingFunction.html
□更に次を予想する。（同表では8桁では5761445）
$ echo "664579" | awk '{print $1*(8.5+0.5)}'
5981211
□次の予想では、1桁前が約3分かかり、10*0.5倍多い比率で素数が見つかるはずなので、
約150分かかるだろうと予測した。
$ echo "10" | awk '{print ($1*3)*($1*0.5) "分 = " ($1*3)*($1*0.5)/60 "時間"}'
150分 = 2.5時間
□現在約74分経過しているので、もしかしたらまた強制終了となるかも知れない。
$ ps -ef | grep clisp | cut -c 40- | awk '{print $1 " " $2}'
01:14:44 /usr/lib/clisp-2.48/base/lisp.run
00:00:00 grep
□が、約89分で終了した。
$ time ./ceron_prime.sh 10000000 `echo "5761455" | awk '{print $1-664579}'` | \
column > 10000000.out
real 89m0.501s
user 88m53.753s
sys 0m20.797s
□上記から合計値=素数の数。特に難しい事をしなくても、
1億までなら正確に求められるということが分かった。
$ wc -w 1*.out
4 1.out
21 10.out
143 100.out
1061 1000.out
8363 10000.out
68906 100000.out
586081 1000000.out
5096876 10000000.out
5761455 合計
□比率についてはほぼ予想通り。
$ wc -w 1*.out | head -8 | awk '{print $1}' | xargs echo -n | \
awk '{print ($1+$2+$3+$4+$5+$6+$7+$8)/($1+$2+$3+$4+$5+$6+$7)}'
8.66933
□一応、素数の数の合計値。
$ wc -w 1*.out | head -8 | awk '{print $1}' | xargs echo -n | \
awk '{print $1+$2+$3+$4+$5+$6+$7+$8}'
5761455
□余談だが経過時間について
$ echo `date '+%H:%M'; ps -ef | grep clisp | grep -v grep | \
awk '{print ":" $5}'` | awk -F\: '{print ($1*60)+$2-($3*60)+$4 "min"}'
99min
--------------------------------------------------------------------------------
■その他、面白そうなもの
参考：P.Graham 著 ANSI Common LISP 練習問題解答
http://www.shido.info/lisp/pacl2.html
