移植
前回のlibm実験の時、OpenBSDにgambit-cをやむなく入れた。だって、gaucheが入らない んだもの。で、折角入れたんだから、一過性のものにしないで何か動かしたい。
そこで考えた。SECDマシンがgauche上で動いていたけど、これをgambit-cに移植して、 C語でコンパイルして単独アプリが出来ないかと。
いきなり、コンパイルして、足りない所を教えて貰おう。
[ob: tcall36]$ gsc -exe vm.scm
*** WARNING -- "list*" is not defined,
*** referenced in: ("/home/sakae/secd/tcall36/vm.c")
*** WARNING -- "list-copy" is not defined,
*** referenced in: ("/home/sakae/secd/tcall36/vm.c")
*** WARNING -- "pop!" is not defined,
*** referenced in: ("/home/sakae/secd/tcall36/vm.c")
*** WARNING -- "push!" is not defined,
*** referenced in: ("/home/sakae/secd/tcall36/vm.c")
ワーニングを出しつつも、vmって名前のアプリが出来上がった。果たして動くか？ 半信半疑ながら実行してみる。
[ob: tcall36]$ ./vm >>> (define aa 1) *** ERROR IN comp -- Operator is not a PROCEDURE (#!unbound 'def 'aa '(stop))
そして、うな垂れる。確信犯なんだけど、コンパイラーは何故エラーにせずにワーニングという 弱気な警告しか出さなかったの？ 多分複数のファイルをコンパイルして、それを合体 する事もあろうに。そういう時は、一つのファイルだけに注目して、関数が定義してない ってエラーには落とせない。だから、せめて警告だけでもって事か。
まあ、いずれにしても、足りないものを補う必要が有る。gaucheのマニュアルを 調べると、list*は、srfi-1のcons*と同じもの。push/popはcommonLispから借用して きたとか書いてある。list-copyはsrfi-1に定義される、便利なやつなんだけど、gambit には無いやつ。
そういう調べの元、srfi-1.scm から、必要な物を頂いてくる。残るは、push/popだけど、これはもう、マクロだろうって 当りをつける。
果たして、gambitでマクロ機能なんて使えるの？ こういう時は、infoするんだった。 マクロ定義が出来ると書いてあった。例でpushが出てた。後はpopを用意すればいいなってんで、 追加した。本当はgemsymで保護が必要だけど、取りあえず無視。下記を追加した。
(define (list* first . rest)
(let recur ((x first) (rest rest))
(if (pair? rest)
(cons x (recur (car rest) (cdr rest)))
x)))
(define (list-copy lis)
(let recur ((lis lis))
(if (pair? lis)
(cons (car lis) (recur (cdr lis)))
lis)))
(define-macro (push! var #!optional val)
`(set! ,var (cons ,val ,var)))
(define-macro (pop! var)
`(let ((head (car ,var)))
(set! ,var (cdr ,var))
head))
これで移植完了かな？
[ob: tcall36]$ time gsc -exe vm.scm
0m9.48s real 0m2.64s user 0m6.48s system
[ob: tcall36]$ ./vm
>>> (define fib (lambda (n) (if (< n 3) n (+ (fib (- n 1)) (fib (- n 2))))))
ldfld00ldc3args2ldg<appselrld00rtnld00ldc1args2ldg-appargs1ldgfibappld00ldc2args2ldg-appargs1ldgfibappargs2ldg+tapprtndeffibstop(time (vm '() '() expr '()))
0 ms real time
0 ms cpu time (0 user, 0 system)
no collections
264 bytes allocated
1 minor fault
9 major faults
fib
>>> (fib 30)
ldc30args1ldgfibappstop(time (vm '() '() expr '()))
8174 ms real time
8170 ms cpu time (5150 user, 3020 system)
1897 collections accounting for 518 ms real time (40 user, 380 system)
2036832352 bytes allocated
289 minor faults
4 major faults
1346269
コンパイラー版だから、スイスイ動くかと思ったら、のっそりした動きだな。
[ob: tcall36]$ gsi vm.scm
>>> (fib 30)
ldc30args1ldgfibappstop(time (vm '() '() expr '()))
146575 ms real time
146520 ms cpu time (122090 user, 24430 system)
13084 collections accounting for 9559 ms real time (1980 user, 7870 system)
10964420600 bytes allocated
225 minor faults
no major faults
1346269
インタープリタは押して知るべしだ。でも、traceを自前で持ってるぞ。
> (trace vm) > (repl) >>> (define hoge (lambda (x) (* x x x))) | > (vm '() '() '(ldf (ld (0 . 0) ld (0 . 0) ld (0 . 0) args 3 ldg * app rtn)... | > (vm '((closure (ld (0 . 0) ld (0 . 0) ld (0 . 0) args 3 ldg * app rtn) ()... | > (vm '(hoge) '() '(stop) '()) | hoge hoge >>> (hoge 5) | > (vm '() '() '(ldc 5 args 1 ldg hoge app stop) '()) | > (vm '(5) '() '(args 1 ldg hoge app stop) '()) | > (vm '((5)) '() '(ldg hoge app stop) '()) | > (vm '((closure (ld (0 . 0) ld (0 . 0) ld (0 . 0) args 3 ldg * app rtn) ()... | > (vm '() '((5)) '(ld (0 . 0) ld (0 . 0) ld (0 . 0) args 3 ldg * app rtn) '... | > (vm '(5) '((5)) '(ld (0 . 0) ld (0 . 0) args 3 ldg * app rtn) '((() () (s... | > (vm '(5 5) '((5)) '(ld (0 . 0) args 3 ldg * app rtn) '((() () (stop)))) | > (vm '(5 5 5) '((5)) '(args 3 ldg * app rtn) '((() () (stop)))) | > (vm '((5 5 5)) '((5)) '(ldg * app rtn) '((() () (stop)))) | > (vm '((primitive #<procedure #2 *>) (5 5 5)) '((5)) '(app rtn) '((() () (... | > (vm '(125) '((5)) '(rtn) '((() () (stop)))) | > (vm '(125) '() '(stop) '()) | 125 125
そしてdebug関連も充実してる。(help-browser "w3m") (help break)とかでも良い。 w3mを指定せずに使いたい時は、/usr/local/bin/gambc-doc を事前に編集しておく。
> ,? ,? : Summary of comma commands ,h | ,(h X) : Help on procedure of last error or procedure/macro named X ,q : Terminate the process ,qt : Terminate the current thread ,t : Jump to toplevel REPL ,d : Jump to enclosing REPL ,c | ,(c X) : Continue the computation with stepping off ,s | ,(s X) : Continue the computation with stepping on (step) ,l | ,(l X) : Continue the computation with stepping on (leap) ,N : Move to specific continuation frame (N>=0) ,N+ | ,N- : Move forward/backward by N continuation frames (N>=0) ,+ | ,- : Like ,1+ and ,1- ,++ | ,-- : Like ,N+ and ,N- with N = nb. of frames at head of backtrace ,y : Display one-line summary of current frame ,i : Display procedure attached to current frame ,b | ,(b X) : Display backtrace of current continuation or X (cont/thread) ,be | ,(be X) : Like ,b and ,(b X) but also display environment ,bed | ,(bed X) : Like ,be and ,(be X) but also display dynamic environment ,e | ,(e X) : Display environment of current frame or X (proc/cont/thread) ,ed | ,(ed X) : Like ,e and ,(e X) but also display dynamic environment ,st | ,(st X) : Display current thread group, or X (thread/thread group) ,(v X) : Start a REPL visiting X (proc/cont/thread)
ああ、emacs用のコードも提供されてるんで、(require 'gambit)しておいて、run-scheme で行ける。
gambit vs. gauche
オイラーの感覚では、gsc版がもっと速いと思っていたんだけど、どうなのよ？ こういう時は、同じ土俵で比べるのがベターです。この俺がルールだとか、 ルールを勝手に変更してしまう某スポーツ界みたいなことは、いけません。
舞台をFreeBSDに移して、gaucheと比較します。 まずは、基準走のgauche君から。
>>> (fib 30) (ldc 30 args 1 ldg fib app stop) ;(time (vm '() '() expr '())) ; real 24.533 ; user 24.352 ; sys 0.078 1346269
gscでコンパイルした、SECDマシン
>>> (fib 30)
ldc30args1ldgfibappstop(time (vm '() '() expr '()))
9717 ms real time
9682 ms cpu time (8411 user, 1271 system)
1944 collections accounting for 881 ms real time (867 user, 587 system)
2036832352 bytes allocated
242 minor faults
no major faults
1346269
gsiのインタープリタ版
>>> (fib 30)
ldc30args1ldgfibappstop(time (vm '() '() expr '()))
127598 ms real time
127005 ms cpu time (119726 user, 7278 system)
13549 collections accounting for 12405 ms real time (11710 user, 3872 system)
10963674080 bytes allocated
179 minor faults
no major faults
1346269
スピードを稼ぐなら当然のごとく、コンパイル版。debugの充実を求めるならgsi。 バランスがよく、何でも有りがgaucheと、予想通りの結果となりました。
cexpl
前回、FreeBSDでオイラーの式を精度よく計算しようとしたら、エラーで刎ねられた。
$ cc t.c -lm
t.c:11:8: warning: implicitly declaring library function 'cexpl' with type
'_Complex long double (_Complex long double)'
z = cexpl(x);
^
t.c:11:8: note: please include the header <complex.h> or explicitly provide a
declaration for 'cexpl'
1 warning generated.
/tmp/t-fd3182.o: In function `main':
t.c:(.text+0x90): undefined reference to `cexpl'
cc: error: linker command failed with exit code 1 (use -v to see invocation)
これが、clangでの結果。曖昧だから、ちゃんと宣言してね。どうかcomplex.hを使って くれってともご宣託。そのくせスパッとエラーに落とされた。
$ gcc48 t.c -lm
t.c: In function 'main':
t.c:11:8: warning: incompatible implicit declaration of built-in function 'cexpl' [enabled by default]
z = cexpl(x);
^
/tmp//ccBpR7vY.o: 関数 `main' 内:
t.c:(.text+0x7a): `cexpl' に対する定義されていない参照です
collect2: error: ld returned 1 exit status
gccでも、ほぼ同様な事を言ってきた。libmがcexpl()を提供してないんだから、当然の 報い。ならば、OpenBSDから、輸血してもらえばいいんでないかい。あれ？ 輸血って 広義の移植なんでしょうか？ ＞医療関係者殿
オイラー的には、血も大事な臓器の一種だと思うから、移植に同意。もしものために 備えて、貯血しておきましょう。人様から血を分けて貰って、同時に肝炎とかHIVも 分けて貰ったんじゃ、本末転倒ですから。
ってな事で、OpenBSDから、分けてもらう。分けて貰った結果は下記。
1 #include <stdio.h>
2 #include <complex.h>
3 #include <math.h>
4
5 long double complex
6 cexpl(long double complex z) {
7 long double complex w;
8 long double r;
9
10 r = expl(creall(z));
11 w = r * cosl(cimagl(z)) + (r * sinl(cimagl(z))) * I;
12 return (w);
13 }
14
15 int main(void){
16 long double _Complex z, x;
17
18 x = I * 2 * acosl(0);
19 z = cexpl(x);
20
21 printf("real=%g\timag=%38.36g\n", creal(z),cimag(z));
22
23 return 0;
24 }
コピペしたんで、本当はクレジットも表記しないとイケンけど許してね。このままだと、 ただのコピペ少年ならぬ、コピペ老人と思われちゃうんで、コードを観賞してみる。 10行目で、実部のべき乗を得る。11行目には、オイラーさんが潜んでいる。cimagl(z)は、 偏角になるんで、オイラー式の一般形となる。
今度は、２人の執刀医による確認でも、文句無く合格。実行もちゃんと出来た。 移植成功で万歳三唱ですかね。
$ cc -g t.c -lm $ gcc48 -g t.c -lm $ ./a.out real=-1 imag=1.2246467991473532071737640294583966e-16
まてまて、結果を注意深く点検するんだ。OpenBSDとかFedoraの結果と比べてミレ。
real=-1 imag=-5.01655761266833226942395050995492469e-20 ;; OpenBSD real=-1 imag=-5.01655761266833226942395050995492469e-20 ;; Fedora real=-1 imag=1.2246467991473532071737640294583966e-16 ;; Debian/ARM
これ、前回の結果から引っ張ってきた。仮数も指数もDebian/ARM 症候群に なってますねぇ。似非、移植が行われたって事？
こんな事もあろうかと、gdb用のプローブを埋め込んでおいたのさ。迷医は、もしも の事も考えておくものです。ブラック・ジャックみたいな天才医はどうなんでしょ。 また、読んでみたいなーー。
で、早速、比較健闘もとえ検討です。医学も正常なものをたくさん集めておいて、対象と 比べる事により進化してますからね。
まずは、正常系のOpenBSD
Breakpoint 1, cexpl (z=0 + 3.1415926535897932385128089594061862 * I)
at /usr/src/lib/libm/src/s_cexpl.c:66
次は異常なFreeBSD
Breakpoint 1, cexpl (z=-1.3222456819757518088907606594413914e+861 + 0 * I)
at t.c:10
計算すべき値がcexplに届いていません。ならば、cexplへの引数を直接指定してみます。 OpenBSDと比較出来るようにcexplはhogeと改名しておきます。
x = 0.0 + M_PI * I; z = hoge(x);
これで、clangでコンパイルしてgdbにかけてみると
Breakpoint 1, hoge (
z=2.4473971540385850596540858107266629e-4942 + <invalid float value> * I)
at t.c:10
10 r = expl(creall(z));
(gdb) up
#1 0x0804879a in main () at t.c:19
19 z = hoge(x);
(gdb) p x
$1 = 0 + 3.1415926535897931159979634685441852 * I
どうも、引数が正しく渡っていないようです。ならば、同じ事をgcc48でやってみると 、正しく渡ってきました。
Breakpoint 1, hoge (z=0 + 3.1415926535897931159979634685441852 * I) at t.c:10 10 r = expl(creall(z));
結局、clangのBugを踏んだって事かな。検証は、糞石レベルまで下りないと出来ない から、これ以上は止めておこう。後は、精度問題。
上で与えてるM_PIでは、精度が取れないようで、OpenBSDでやると、e-16の精度に なる。そこで、x = 0.0 + (2 * acosl(0)) * I; を復活させる。
で、じっとcexplを眺めてみると、虚部に関する精度の肝は、sinlに委ねられているって 予想される。よって、その部分を取り出して検証出来るようにしておく。
long double complex
hoge(long double complex z) {
long double complex w;
long double ci, sv, r;
ci = cimagl(z);
sv = sinl(ci);
:
z=0 + 3.1415926535897931159979634685441852 * I ;; FreeBSD
z=0 + 3.1415926535897932385128089594061862 * I ;; OpenBSD
^^^^^^^^^^^^^^^^^^^
z=0 + 3.1415926535897931159979634685441852 * I ;; 0.0 + M_PI * I
(gdb) p ci
$1 = 3.1415926535897931159979634685441852 ;; FreeBSD
$1 = 3.1415926535897932385128089594061862 ;; OpenBSD
^^^^^^^^^^^^^^^^^^^
(gdb) p sv
$2 = 1.2246467991473531772014792699320386e-16 ;; FreeBSD
$2 = -5.0165576126683320234517576003912636e-20 ;; OpenBSD
どうも、OpenBSDでは、acoslの精度が悪いようだ。そして、そのいかがわしい値を 使って計算した結果が、たまたま誤差が少なくなったものだから、安心しちゃった んだな。
まとめると、sinl(π)を正確に求めましょ。言い換えるとsin(180°)って事だから、 期待値は、中学生でも知ってる答え ゼロ。そのπをacosから導出したものだから、 正しいπにならなかった。誤差が相殺して、たまたま精度よく答えが得られたと 錯覚した。fedoraもOpenBSDと同じ結果だったので、多数決も味方したって事だな。
いやー、オイラーさんは色々な事を勉強させてくれるな。
clangの弁明
上で、clangを使ってコンパイルした時、cexplに引数が正しく渡っていないってサジを 投げてしまった。でも、改めて見直すと、どうやらそうとも言えないと思える。
$ cc -g -O0 t.c -lm
$ gdb -q a.out
Reading symbols from a.out...done.
(gdb) b hoge
Breakpoint 1 at 0x8048692: file t.c, line 10.
(gdb) run
Starting program: /usr/home/sakae/z/a.out
Breakpoint 1, hoge (z=<invalid float value> + <invalid float value> * I)
at t.c:10
10 ci = cimagl(z);
(gdb) n
11 sv = sinl(ci);
(gdb)
13 r = expl(creall(z));
(gdb)
14 w = r * cosl(cimagl(z)) + (r * sinl(cimagl(z))) * I;
(gdb) p z
$1 = <optimized out>
(gdb) p ci
$2 = 3.1415926535897931159979634685441852
(gdb) p sv
$3 = 1.2246467991473531772014792699320386e-16
(gdb) p r
$4 = 1
(gdb) c
Continuing.
real=-1 imag=1.2246467991473532071737640294583966e-16
hoge(cexpl)に渡ってくる引数は、clang流に解釈すれば定数。だから、関数内ではオプチマイズ されててそのままでは参照出来んと言ってるんだな。gdbでbreakした時は、gdbの 解釈に齟齬があって、invalidって言ってるのだろう。 cimaglではきちんと値が取り出せている。また、creallの値は直接には確認して いないけど、計算結果から容易に0と推測される。
ちょっと、objdumpして、該当箇所を覗き見してみる。
x = 0.0 + M_PI * I; 8048763: db bc 24 98 00 00 00 fstpt 0x98(%esp) 804876a: dd 05 d0 88 04 08 fldl 0x80488d0 8048770: d9 c0 fld %st(0) 8048772: db bc 24 a4 00 00 00 fstpt 0xa4(%esp) z = hoge(x); 8048779: db ac 24 98 00 00 00 fldt 0x98(%esp) 8048780: db 7c 24 68 fstpt 0x68(%esp) 8048784: db 7c 24 74 fstpt 0x74(%esp) 8048788: 8b 4c 24 7c mov 0x7c(%esp),%ecx :
インテルのマニュアルから動作を確認しようとしたけど微妙にニューモニックが異なる。 でも、定数を直接x87に転送してるっぽい。動作の確認は面倒なので、gdbに任せちゃえ。
Breakpoint 1, main () at t.c:21
21 x = 0.0 + M_PI * I;
(gdb) n
22 z = hoge(x);
(gdb) info float
=>R7: Valid 0x4000c90fdaa22168c000 +3.141592653589793116
R6: Empty 0x4000c90fdaa22168c000
R5: Empty 0x00000000000000000000
R4: Empty 0x00000000000000000000
R3: Empty 0x00000000000000000000
R2: Empty 0x00000000000000000000
R1: Empty 0x00000000000000000000
R0: Empty 0x00000000000000000000
Status Word: 0x3800
TOP: 7
Control Word: 0x127f IM DM ZM OM UM PM
PC: Double Precision (53-bits)
RC: Round to nearest
Tag Word: 0x3fff
Instruction Pointer: 0x33:0x08048772
Operand Pointer: 0x3b:0xbfbfe804
Opcode: 0xdbbc
(gdb) display/i $pc
1: x/i $pc
=> 0x8048779 <main+57>: fldt 0x98(%esp)
思った通り、0x08048772番地の命令(hogeへの突入直前)を実行した後、パイの値が FPUに転送されている。これじゃ、gdbが右往左往したってしょうがないか。虚部の パイがFPUの載るのは良いとして実部の0はどうするんだろう？ そこがちと疑問。
hogeの中をちょっと変更してから、コードを覗いてみる。
hoge(long double complex z) {
8048680: 55 push %ebp
8048681: 89 e5 mov %esp,%ebp
8048683: 81 e4 f8 ff ff ff and $0xfffffff8,%esp
8048689: 81 ec 80 00 00 00 sub $0x80,%esp
804868f: 8d 45 0c lea 0xc(%ebp),%eax
8048692: 8b 4d 08 mov 0x8(%ebp),%ecx
long double complex w;
long double ci, cr, sv, r;
cr = creall(z);
8048695: db 28 fldt (%eax)
8048697: db 7c 24 50 fstpt 0x50(%esp)
ci = cimagl(z);
804869b: db 68 0c fldt 0xc(%eax)
804869e: d9 c0 fld %st(0)
80486a0: db 7c 24 5c fstpt 0x5c(%esp)
どうやらzのデータは、ポインターで渡ってきているな。そして実部も虚部も、 一度FPUに載せて、それをfstptで、所定のローカルエリアに転送してる。 こういうのも有りなんだろうけど、そんな作法はGNUの世界では一般的じゃ ないので面倒見れんと言う事か。あー、すっきりした。
FreeBSDのデフォルトccがclangになって、皆さんdebugには苦労していないんですかね？ オイラーみたいな軟弱者のために、clangとタイアップしたlldbがリリース されないかな？
後は、折角オイラーさんが、libmの中に潜んでおられるので、よくそのお姿を 拝んでおけ。