bc/dc
ニュートン誌より
中坊？向けの総合科学雑誌として、もう40年もニュートンが発行されてる。 数学記事も度々顔を出す。人気は、素数、パイ、虚数が定番。加えて、三角関数、指数・対数、微分・積分あたりだろうか。
6月号は、数学教養教室って名うってたから、一般人向けとも考えられるな。で、お題は、指数・対数。今回は、どんな切りぐち？
昔、功績をあげた小僧が、殿様からほうびをつかわそう、と言われた。それに対して、今日は米粒１ヶ、明日はその倍に２こ、明後日は、その倍の４個、、これを30日続けてください。なんじゃ、そんな亊でいいのか。。こういう古典過ぎるのは流行らない。
今回は、フルマラソン目指してのトレーニング。最初から42Kmなんて無理に決ってる。最初は1Kmからのスタート。翌日は、それより1%増し、、、って具合に、無理なく1年続けるのじゃ。 で、1年後、何Kmを走る練習になってるかな？
福利計算ともいえるな。1.01を365回かけ算するといいよ。馬鹿みたいに電卓を叩きますか？
a^m X a^n = a^(m+n)
あたりを参考に、電卓を叩く回数を減らす工夫をしてみる？ オイラーは、文明の利器を利用して、答一発ですよ。
gosh> (expt 1.01 365) 37.78343433288728 gosh> (expt 1.01 377) 42.57531953064501
１年だと、わずかに足りない。もうちょっと頑張れば、なんとかなる。
昔の忍者は、跳躍の練習の為、麻の種を植えたそう。芽が出たら、それをひたすら飛び越える、あけてもくれても飛び越える。伸てきても飛び越える。我慢がきつくなったら、麻でハイになって、超人力を発揮する、らしいですよ。
dc
前回の最後にやった、謎のコマンド列を紐解いてみたい。一応、実行例を再掲。10までの階乗を求めている。
[la1+dsa*pla10>y]sy 0sa1 lyx 1 2 6 24 : 3638800
OpenBSDの説明書は簡潔過ぎたので、リナのdc(1)を参考にした。対象毎に整理されてるので、理解が進むだろう。下記は、関係ありそうなものを見繕っておいた。
Registers
sr Pop the value off the top of the stack and store it into regis‐
ter r.
lr Copy the value in register r and push it onto the stack. The
value 0 is retrieved if the register is uninitialized. This
does not alter the contents of r.
Stack Control
c Clears the stack, rendering it empty.
d Duplicates the value on the top of the stack, pushing another
copy of it. Thus, ``4d*p'' computes 4 squared and prints it.
Strings
a The top-of-stack is popped. If it was a number, then the low-
order byte of this number is converted into a string and pushed
onto the stack. Otherwise the top-of-stack was a string, and
the first character of that string is pushed back.
x Pops a value off the stack and executes it as a macro. Normally
it should be a string; if it is a number, it is simply pushed
back onto the stack. For example, [1p]x executes the macro 1p
which pushes 1 on the stack and prints 1 on a separate line.
>r Pops two values off the stack and compares them assuming they
are numbers, executing the contents of register r as a macro if
the original top-of-stack is greater. Thus, 1 2>a will invoke
register a's contents and 2 1>a will not.
Printing Commands
p Prints the value on the top of the stack, without altering the
stack. A newline is printed after the value.
f Prints the entire contents of the stack without altering any‐
thing. This is a good command to use if you are lost or want to
figure out what the effect of some command has been.
コマンドは、基本が一文字だ。中には、例外があって、それはレジスターとスタック間のデータ転送とかになる。レジスター名前は、通常一文字(-xをつけてdcを起動すると拡張出来る)。
レジスターって言うけど、様は変数名だ。それから、文字列は、[srring] のように定義する亊になっている。
暗号解読
さあ、暗号解読だ。一行目は、文字列をスタックに配置してから、それをyレジスターに保存。 bc言語風にかくと、define y(..){ .. } だな。
関数だから、多分引数が必要だろう。引数はどうやって渡す？ 根がスタック言語(forthの劣化版もとえ、簡略版)だから、スタックに積んで渡す。あるいは、レジスタターと称する変数に入れて渡すの二種類が考えられるな。
そんな亊を思いながら、先に3行目を見る。lyx って、関数の実行そのものじゃん。先に関数をyレジに保存した。それを呼び戻して(stackに積んで)それをxで実行。
となると、2行目は、引数の準備だな。面倒だから実行してみるか。
0sa1 f 1 la f 0 1
0saって、a=0 だな。次に数値の1がスタックに詰れる。この時点でfを実行し、スタックの全ダンプをした。勿論、結果は1だ。念の為、laで、aレジスターを呼び戻してみると、当然0が詰れている。
以上の結果から、引数の一つはaレジスターに、もう一つはスタックに詰れている亊が判明した。また、上記を実行した後、スタックを確認したら、最後の結果が残っていた。
この事実を考慮すると、 スタックは結果を保持するんだろう。初期値が１っての頷ける。それからaレジスターは、何回実行したかのカウンターの役を仰せ使っているのだろう。
以下で、関数を表す文字列を、関数っぽくデコードしてみる。一応スタックのトップを大域変数と考えて、関数を呼び出した時の値が入ってる変数名をTって亊にしとく。Uはもう一つのスタックの値とする。
２行目と３行目は、次のようになるな。
T = 1 a = 0 y(a)
define y(a) {
T
T U = 1 a + ;; la1+
T U U ;; d
a = U ;; sa ( a = a + 1)
T U ;; saで、１ヶ消費された
T = T * U ;; *
T p ;; print
T a 10 > ;; la10> aレジと10の比較。
T call y(a) ;; a < 10 なら、再度、関数yを呼出
}
中々、スタック言語を変数で表すのは難しい。が、２度目に呼び出される時も、最初の呼出形態と同じになってる。但し、スタックにはかけ算した結果が、引数は１増えた状態での呼出になってる。
try newton
gaucheは普通の人は持っていないだろう。でも、dcぐらいは入ってると仮定して、それで計算させてみる。
上の暗号を解いたおかげで、すらすらと書けたよ。 p.dc
[la * lc 1- d sc 0<y]sy 20k 1.01sa 377sc 1 lyx p
実行結果は、377日の成果ね。
vbox$ dc p.dc 42.57531953064487204540
関数を定義、呼び出す準備、呼び出して、結果がスタック残るのでそれを表示。
aレジスタに明日の進歩度合い、cレジスタに努力日数、１は初期値。ああ、20kは、少数点以下の桁数。
関数内は、分かり易くなるように、適当にスペースを入れている。スタックの内容とaレジを積算、cレジを1減算、結果を複製、cレジに退避、最後は、日数が残っていたら、再度関数を呼出。関数から抜けた時、スタックに努力結果が残っている。こんな塩梅ね。
[la * lc 1- d sc 0<p]sp 1.01sa [How many days ]P ? sc 20k 1 lpx p
? を使うと、stdinから入力出来るようになる。
vbox$ dc p.dc How many days 365 37.78343433288715886170
ちなみに、上で出て来た377日ってのは、どうやって出てきた。ニュートン誌からのコピー？ なんて言うと中坊って馬鹿にされるぞ。
1.01^n = 42.195 のnを求めるのだ。 l(42.195)/l(1.01) 376.09822083638752008469
答を切り上げて377日ね。任意の底の対数を求める変換式を使ってみた(下記のcontrib参照)。
yydebug for bc
OpenBSDのソースを見ていたら、やたYYDEBUGってのが目についたので、何かご利益あるだろうってんで、探りを入れてみた。yaccの結果ファイル、 bc.c の方ね。
#ifndef YYDEBUG
#define YYDEBUG 1
#endif
:
yyparse(void)
{
int yym, yyn, yystate;
#if YYDEBUG
const char *yys;
if ((yys = getenv("YYDEBUG")))
{
yyn = *yys;
if (yyn >= '0' && yyn <= '9')
yydebug = yyn - '0';
}
#endif /* YYDEBUG */
でソースコードにdebugコードを注入し、export YYDEBUG=1 で、enableすると言う、慎重さを潜り抜けると、
vbox$ ./bc -c yydebug: state 0, reducing by rule 1 (program :) yydebug: after reduction, shifting from state 0 to state 1 12345 + 6789 yydebug: state 1, reading 268 (NUMBER) yydebug: state 1, shifting to state 7 yydebug: state 7, reducing by rule 72 (expression : NUMBER) yydebug: after reduction, shifting from state 1 to state 28 yydebug: state 28, reading 296 (PLUS) yydebug: state 28, shifting to state 73 yydebug: state 73, reading 268 (NUMBER) yydebug: state 73, shifting to state 7 yydebug: state 7, reducing by rule 72 (expression : NUMBER) yydebug: after reduction, shifting from state 73 to state 111 yydebug: state 111, reading 266 (NEWLINE) yydebug: state 111, reducing by rule 76 (expression : expression PLUS expression) yydebug: after reduction, shifting from state 1 to state 28 yydebug: state 28, reducing by rule 19 (statement : expression) yydebug: after reduction, shifting from state 1 to state 33 yydebug: state 33, reducing by rule 8 (semicolon_list : statement) yydebug: after reduction, shifting from state 1 to state 32 yydebug: state 32, shifting to state 83 yydebug: state 83, reducing by rule 3 (input_item : semicolon_list NEWLINE) 12345 6789+ps. yydebug: after reduction, shifting from state 1 to state 31 yydebug: state 31, reducing by rule 2 (program : program input_item) yydebug: after reduction, shifting from state 0 to state 1
こんな楽しい情報が得られる。なお、bcには全く関係無いと思われる、/usr/lib/libcrypto.so.49.0なんてのをリンクしてた。何でかなと思ったら、こいつが提供してる多倍長ルーチンを使いたかったからと判明。ほとんどdcのコードそのものを使ってた。
bc of contrib
FreeBSDのbcが寄贈品と言う亊は、寄贈元が有る訳で、調べてみると、 gavin/bc だった。Windowsも全々大丈夫だよ、イェーいって所かな。でも、Windowsユーザーに黒い窓は受入られるのかな？ 最近、LinuxにGUIは甘えなのか こんな記事が出てるからなあ。リナな人も、黒い窓怖いって人が多いだろうに。事実、リナでもbc/dcはデフォでは入っていません。
で、つらつらソースを眺めていたら楽しい機能を発見。ニュートン誌の例題を、一発で解いてくれる関数が提供されてた。
[sakae@fb ~]$ bc -l p(1.01, 365) 37.78343433288715887761 p(1.01, 377) 42.57531953064487206341
これ、POSIXに定義されてるbcを拡張してるのね。ソースは、 gen/lib2.bc
define p(x,y){
auto a
a=y$
if(y==a)return (x^a)@scale
return e(y*l(x))
}
大事なのは、最後のeとlを使った計算。lは、ナチュラル・ロガリズム。eを底にした対数だ。
ニュートン誌には、このナチュラル・ロガリズムを使って、任意の数を底にしたログを求める公式が紹介されてる。log a b = l(b) / l(a) てやつ。aは任意な底、bは求めたい数。(注: lはbc上の関数名、数学では、ln と表記される亊が多い)
これを使えば、簡単にlog 10 x が計算出来る。log 20 = l(20) / l(10) って具合。
上のpって関数は、このログの計算と対になる関数だな。こんなの初めて知ったぞ。
リナでGNUに毒されていないbcをやってみる。
昔、あのペンギン好きな人が趣味でカーネルを作っていた。、リチャード・ストールマンがみはてぬ夢のフリーなOSを作っていた。unixは金まみれ、権利まみれで自由にならないって理由からだ。
所がさっぱり、このフリーなOSが完成しない。ええい、面倒だ。フリーなOSのカーネル部分をペンギン製に置き換えちゃえ。そしてOSの名前はペンギン好きさんの名前を頂いてリナックスにしよう、ですからね。だから、カーネル以外は全部GNU謹製なのさ。
sakae@deb:/tmp/bc$ ls -l bin/*
-rwxr-xr-x 1 sakae sakae 181928 Jun 11 06:44 bin/bc*
lrwxrwxrwx 1 sakae sakae 2 Jun 11 06:44 bin/dc -> bc*
sakae@deb:/tmp/bc$ ldd bin/bc
linux-gate.so.1 (0xb7f91000)
libc.so.6 => /lib/i386-linux-gnu/libc.so.6 (0xb7d5d000)
/lib/ld-linux.so.2 (0xb7f93000)
変なライブラリィーを全く使わない、移植のための異色なアプリです。自前で多倍長ルーチンも用意してるのね。
history of bc/dc
ふと思いたって、bc/dc比べをしてみた。
ob$ ls -l /usr/bin/{bc,dc}
-r-xr-xr-x 1 root bin 100848 Apr 12 08:46 /usr/bin/bc*
-r-xr-xr-x 1 root bin 45440 Apr 12 08:46 /usr/bin/dc*
サイズが倍以上違うな。bcはパース部分が入っているからなあ。dcは、早い話RPNな計算器だ。 以前やったmrubyがbcに対応。そこから実行機能だけ取出したのがmruby/cだった。dcに相当するんだな。bc側からみたら、仮想マシーンって亊になる。
そうなると、bcが先に出来た？ それともdcが先、それとも同時開発？ 一応manをあたってみる。
HISTORY
The bc command first appeared in Version 6 AT&T UNIX. A complete rewrite
of the bc command first appeared in OpenBSD 3.5.
HISTORY
The dc command appeared in Version 1 AT&T UNIX. A complete rewrite of
the dc command using the BN_new(3) big number routines first appeared in
OpenBSD 3.5.
答はdcが先でした。
所で、version 6 AT&AT UNIXっていつ頃の亊？ Unixの歴史 に年表が出てた。1975年とは時代がかっているな。 久し振りに昔に戻ってみるか。ええと、今迄何度かやっているな。 ニュートンみたいに毎年じゃないから許せ。
ソースは、usr/source/s1 の中。chdir で移動する。
# ls -l bc* dc* -rw-r--r-- 1 bin 11348 Jan 1 1970 bc.y -rw-r--r-- 1 bin 1773 Jan 1 1970 bcd.c -rw-r--r-- 1 bin 29076 Jan 1 1970 dc1.s -rw-r--r-- 1 bin 6009 Jan 1 1970 dc2.s -rw-r--r-- 1 bin 5885 Jan 1 1970 dc3.s -rw-r--r-- 1 bin 5396 Jan 1 1970 dc4.s -rw-r--r-- 1 bin 7432 Jan 1 1970 dc5.s -rw-r--r-- 1 bin 3332 Jan 1 1970 dcheck.c
bcの元になるbc.yでの供給ですか。それからdcは、pdp11のアセンブラーによる努力の結晶です。この頃のOSはhead,tail,wcも無い超越したもの。こんな環境で、よく開発出来たものだ。頭が下がりますよ。
ちょいと、bc.yを覗いてみる。
# cat bc.y
%right '='
%left '+' '-'
%left '*' '/' '%'
%right '^'
%left UMINUS
%term LETTER DIGIT SQRT _IF FFF EQ
%term _WHILE _FOR NE LE GE INCR DECR
%term _RETURN _BREAK _DEFINE BASE OBASE SCALE
%term EQPL EQMI EQMUL EQDIV EQREM EQEXP
%term _AUTO DOT
%term QSTR
stat : e
={ bundle( $1, "ps." ); }
|
={ bundle( "" ); }
| QSTR
={ bundle("[",$1,"]P");}
:
e : e '+' e
= bundle( $1, $3, "+" );
| e '-' e
= bundle( $1, $3, "-" );
| '-' e %prec UMINUS
= bundle( " 0", $2, "-" );
これを見ると、bc上で、1+2 ってやったのが、1 2 + ps. って、翻訳される亊が分る。で、翻訳結果はdcが受け取るはず。 bc.yのmainを見ると、
pipe(p);
if (fork()==0) {
close(1);
dup(p[1]);
close(p[0]);
close(p[1]);
yyinit(argc, argv);
yyparse();
exit();
}
close(0);
dup(p[0]);
close(p[0]);
close(p[1]);
execl("/bin/dc", "dc", "-", 0);
forkして、子供側はbcを実行。親側はdcに挿げ替えて、パイプで翻訳結果をdcに送るとな。見事な連携であります。
また、/usr/lib/lib.b に拡張ファイルが置いてあった。
# bc -l e(1) 2.71828182845904523536
勿論、ニュートンの例題も、問題なく実行できたよ。
parseって、超昔から有る有名なやつなんですねぇ。以上、軽く歴史散策してみました。
ついでに各dirにrunってのが配置されてて、これがコンパイルの手順を表している。該当部分と他の例
yacc bc.y cc -s -O y.tab.c -ly cmp a.out /usr/bin/bc cp a.out /usr/bin/bc rm y.tab.c cc -s -O cal.c cmp a.out /usr/bin/cal cp a.out /usr/bin/cal cc -s -n -O cc.c cmp a.out /bin/cc cp a.out /bin/cc
次は、思いっきり現代へ飛躍します。
compiler
と、言う訳で、新しいコンパイラーは、どうやってるか、調べてみる。
これ素晴しい。yacc無しで自前で解析方法を構築している。また、随所にこぼれ話が挿入されてて楽しいぞ。 例えば、今のC言語の不合理さは歴史の所産とか、昔のコンパイラーがどうなってたも言及されてる。上で昔にワープしたのが、リアルタイムで観測出来たよ。こういうのが、本当の楽しみなんだね。