dc (2)

以前パソコンを買いたいと相談を受けた人から、新築祝いの報告を受けた。

2年使用の中古品に、Windows10、256GのSDD、クワッドの石なやつが、 52000円とか。

21インチのモニターが17000円、キーボードが1000円。無線LANが親子セットで 7000円。マウスだけは古いのを使用したって言ってた。

本体のイヤホンジャックが接触不良でいらいらしたので、1500円のUSBオーディオを用意したとか。マイクもついてて、これってスカイプするための汎用品と想像するけど、どうよ。USBオーディオもピンキリみたいで、流行のハイレゾ対応とかがパソコンよりも高価だったりするみたい。需要の有る所に供給が有るんだなあ。

数秒で立ち上がるし、無線の親が２Fに有って、子が１Fのリビングに置いても、さくさくと動画が見られると喜んでいた。

パソコン本体はつてをつたって中古を使ったらしいけど、気をきかせて firefoxやらChromeを入れてくれていたので、何もせぬとも、使えるようになったとな。後は、家計簿を付けるのに、オープンオフィスを入れると言っていた。

それから、ワードプレスとか言うのを入れて、ホームページをやりたいとか。なんじゃそれ？

Wordpressとは？ワードプレスのインストール方法・使い方と、オススメの無料テーマ（テンプレート）・プラグイン紹介 phpとMySQLをセットにしたブログシステムか。なんか、恐い人の餌食にされそうですよ。先にこちらを見て、知識を蓄えてからね。

Webアプリの脆弱性調査ツール、脆弱性学習ソフトなど

Linuxのdc

前回、Debian（親戚のウブでも同じ）に、dcを入れたので、軽く見ておく。一番の注目点は、大きな数値をどう扱っているかだ。

FAQに、しっかり言い訳してあった。

1) Why does BC have its own arbitrary precision number routines
   (found in lib/number.c) rather than using GMP?

GMP has "integers" (no digits after a decimal), "rational numbers"
(stored as 2 integers) and "floats".  None of these will correctly
represent a POSIX BC number.  Floats are the closest, but will not
behave correctly for many computations.  For example, BC numbers have
a "scale" that represent the number of digits to represent after the
decimal point.  The multiplying two of these numbers requires one to
calculate an exact number of digits after the decimal point regardless
of the number of digits in the integer part.  GMP floats have a
"fixed, but arbitrary" mantissa and so multiplying two floats will end
up dropping digits BC must calculate.

GMPは帯に短したすきに長しだから、BC用の無限精度ルーチンを書いたよ。車輪の大発明をやっちまっただつう事みたいだな。これはunix精神を逸脱してると思うんだけど、どうだろう？余り勝手な事を書くと、どこからかマサカリが飛んできそうなので、これぐらいにしておく。

sakae@ub:~/testme/bc-1.06.95$ sh debian/dc_factor 12345678901
857
14405693

[ob: ~]$ cat dc_factor
#!/bin/sh

dc -e "$1[p]s2[lip/dli%0=1dvsr]s12sid2%0=13sidvsr[dli%0=1lrli2+dsi!>.]ds.xd1<2"

これ、Debianのおまけ。このコードが読めれば、dcの免許皆伝、かな？それじゃ、dcの中身をちょっと、追ってみる。

(gdb) p/c c
$4 = 53 '5'
(gdb) bt
#0  dc_func (c=53, peekc=100, negcmp=0) at eval.c:175
#1  0x0804a441 in dc_evalstr (string=0xbffff594) at eval.c:572
#2  0x08049271 in main (argc=3, argv=0xbffff654) at dc.c:308

これ、引数に、-e '5d*p' を与えた所。

=>  switch (c){
    case '_': case '.':
    case '0': case '1': case '2': case '3':
    case '4': case '5': case '6': case '7':
    case '8': case '9': case 'A': case 'B':
    case 'C': case 'D': case 'E': case 'F':
        return DC_INT;
    case ' ':
    case '\t':
    case '\n':
        /* standard command separators */
        break;

    case '+':   /* add top two stack elements */
        dc_binop(dc_add, dc_scale);
        break;
    case '-':   /* subtract top two stack elements */
        dc_binop(dc_sub, dc_scale);
        break;

こんな感じで、コマンド解析が行われている。あちらこちらに飛んでて、どうもすっきりしない。すっきりしない、もう一つの理由は、ソースのあちこちで、#ifdefが出て来る事。これに遭遇すると、追いかけていて、そわそわしてしまう。まあ、多システムで動くようにすると、こうなってしまうのは、必然だけど。。。

Fedoraでは、どうよ

リナのもう一つの雄である、Fedora系では、ソースのお取り寄せと解析はどうやればいいの？軽く調べておくかな。

nginxをパッケージとしてコンパイルしていろいろ

どうも、Fedoraの方は面倒くさそうだなあ。

最初にパッケージを作る為の場所を定義。dcを取ってこようとしたけど無いと言われたのでbcを取ってくる。

[sakae@fedora ~]$ echo %_topdir ~/rpm > .rpmmacros
[sakae@fedora testme]$ dnf download --source bc
enabling fedora-source repository
enabling updates-source repository
bc-1.06.95-15.fc23.src.rpm                      847 kB/s | 306 kB     00:00
[sakae@fedora testme]$ ls
bc-1.06.95-15.fc23.src.rpm

そして、展開。すると、上で指定した所に入る。

[sakae@fedora testme]$ rpm -ivh bc-1.06.95-15.fc23.src.rpm
更新中 / インストール中...
   1:bc-1.06.95-15.fc23               警告: ユーザー mockbuild は存在しません - root を使用します
  :
警告: グループ mockbuild は存在しません - root を使用します
################################# [100%]

こんな具合になった。

[sakae@fedora ~]$ tree rpm
rpm
├── SOURCES
│ﾂﾂ ├── bc-1.06-dc_ibase.patch
│ﾂﾂ ├── bc-1.06.95-doc.patch
│ﾂﾂ ├── bc-1.06.95-matlib.patch
│ﾂﾂ ├── bc-1.06.95-memleak.patch
│ﾂﾂ ├── bc-1.06.95-sigintmasking.patch
│ﾂﾂ └── bc-1.06.95.tar.bz2
└── SPECS
    └── bc.spec

rpm関係者を入れ、依存者を処理するための者を入れようとするとエラー。

[sakae@fedora ~]$ dnf -y install rpm-build
[sakae@fedora ~]$ sudo dnf buiddep --spec bc.spec
コマンド「buiddep」が見つかりません。「/bin/dnf --help」を実行してください。
DNF プラグインコマンドが利用できません。"dnf install 'dnf-command(buiddep)'" を 実行してください。
[sakae@fedora ~]$ sudo dnf install 'dnf-command(buiddep)'
パッケージ dnf-command(buiddep) は利用できません。
エラー: Unable to find a match.

エラーを無視して、スペックファイルの有る所でコンパイル

[sakae@fedora SPECS]$ rpmbuild -bb bc.spec

こういう場所に目指すものが出来上がっていた。

/home/sakae/rpm/BUILD/bc-1.06.95/dc

ストリップされていないので、gdbにそのままかけられる。 debian系の方が、手軽な感じがするな。

dc on unix-v6

マニュアルに、dcはunix-v6から表れたと書いてあったので、確かめてみる。以前、unix-v6のソース嫁に凝った時に作っておいたものが有るんで、 OpenBSDで確認。

[ob: pdp11]$ ./boot

PDP-11 simulator V3.9-0
./boot> #!/usr/local/bin/simh-pdp11
Unknown command
Disabling XQ
@unix

login: root
# dc
12d*p
144
11111111111111111111111111111111111111111111111d*p
1234567901234567901234567901234567901234567901209876543209876543209876\
54320987654320987654321
q

うん、確かに存在してる。売りの無限精度も実現してるし、いにしえの人は偉かった！で、どんな構成になってたか？

# wc dc[1-5]*
   2542    5078 dc1.s
    457    1070 dc2.s
    460     998 dc3.s
    408     983 dc4.s
    540    1316 dc5.s
   4407    9445 total
# cat dc2.s
   :
.text
/
sqrt:
        mov     r4,-(sp)
        mov     r3,-(sp)
        mov     r2,-(sp)
        mov     r0,-(sp)
/
/       check for zero or negative
/
        mov     w(r1),r2
        sub     a(r1),r2
/
/       look at the top one or two digits
/
        mov     r1,r3
        jsr     pc,fsfile
        jsr     pc,backspace
        mov     r0,r4
        bit     $1,r2
        bne     2f
        mov     r4,r1
        mul     $100.,r1
        mov     r1,r4
        mov     r3,r1
        jsr     pc,backspace
        add     r0,r4

５本のアセンブラーソースになってた。適当に開いてみたら、pdp-11語で実現した平方根のルーチンだったぞ。平方根の記号にはvが割り当てられている。なんとなく、字形がルート記号に似てるという発想なんだろうね。

2vp
1
5k 2vp
1.41421
30k 2vp
1.414213562373095048801688724209

kに前置する数字で、少数点以下の桁数を自由に設定出来る。初期値は、0になってる。

うんと古いOpenBSDでは、どう書かれていたんだろう？ CVS log for src/usr.bin/dc/bcode.cを見る限りでは、OpenBSD3.5から存在する事になってる。まあ、bn(3)を使って書き直したって言ってるんだから、当たり前だけどね。

FreeBSD 6.4で調べてみたら、bc/dcは、contribに置いてあった。GNUから拝借していたんだな。FreeBSD 10.2では、OpenBSDからの移植版が正しい位置の /usr/src/usr.bin/dcに鎮座してた。BSD三兄弟の麗しき友情。

root

平方根はどうやって求める？えっちらおっちらとpdp-11のコードを読むか？まさかね。OpenBSDのdcで読んでみる。

bsqrt(void)
{
   :
        n = pop_number();
   :
                scale = max(bmachine.scale, n->scale);
                normalize(n, 2*scale);
                x = BN_dup(n->number);
                bn_checkp(x);
                bn_check(BN_rshift(x, x, BN_num_bits(x)/2));
                y = BN_new();
                bn_checkp(y);
                ctx = BN_CTX_new();
                bn_checkp(ctx);
                for (;;) {
                        bn_checkp(BN_copy(y, x));
                        bn_check(BN_div(x, NULL, n->number, x, ctx));
                        bn_check(BN_add(x, x, y));
                        bn_check(BN_rshift1(x, x));
                        if (bsqrt_stop(x, y, &onecount))
                                break;
                }
                r = bmalloc(sizeof(*r));
                r->scale = scale;
                r->number = y;
                BN_free(x);
                BN_CTX_free(ctx);
                push_number(r);

本質部分は上記だ。と言っても、直ぐに理解出来たわけではない。平方根のアルゴリズムとか高速sqrt()関数を参考にさせてもらったぞ。

BN_divの引数が沢山有ってくらくらするので、manから抜粋しとく

     int
     BN_div(BIGNUM *dv, BIGNUM *rem, const BIGNUM *a, const BIGNUM *d,
         BN_CTX *ctx);

     BN_div() divides a by d and places the result in dv and the remainder in
     rem (dv=a/d, rem=a%d).  Either of dv and rem may be NULL, in which case
     the respective value is not returned.  The result is rounded towards
     zero; thus if a is negative, the remainder will be zero or negative.  For
     division by powers of 2, use BN_rshift(3).

GNUのroot

GNUの人(達)は、どうやって平方根を求めているのかな？参考に見ておく。場所は、lib/number.c/bc_sqrt

  /* Calculate the initial guess. */
  if (cmp_res < 0)
    {
      /* The number is between 0 and 1.  Guess should start at 1. */
      guess = bc_copy_num (_one_);
      cscale = (*num)->n_scale;
    }
  else
    {
      /* The number is greater than 1.  Guess should start at 10^(exp/2). */
      :
  /* Find the square root using Newton's algorithm. */
  done = FALSE;
  while (!done)
    {
      bc_free_num (&guess1);
      guess1 = bc_copy_num (guess);
      bc_divide (*num, guess, &guess, cscale);
      bc_add (guess, guess1, &guess, 0);
      bc_multiply (guess, point5, &guess, cscale);
      bc_sub (guess, guess1, &diff, cscale+1);
      if (bc_is_near_zero (diff, cscale))
        {
          if (cscale < rscale+1)
            cscale = MIN (cscale*3, rscale+1);
          else
            done = TRUE;
        }
    }

ニュートン法を使うとな。りんごのニュートンさん？

SICP

そう言えば、SICPにも出てきたな。計算機プログラムの構造と解釈第二版

[ob: ~]$ gsi
Gambit v4.8.4

> (define (sqrt-iter guess x)
    (if (good-enough? guess x)
      guess
      (sqrt-iter (improve guess x)
                 x)))
> (define (improve guess x)
    (average guess (/ x guess)))
> (define (average x y)
    (/ (+ x y) 2))
> (define (good-enough? guess x)
    (< (abs (- (square guess) x)) 0.00001))
> (define (mys x)
    (sqrt-iter 1.0 x))

定義したものとデフォルトのものを比べてみる。

> (mys 2)
1.4142156862745097
> (sqrt 2)
1.4142135623730951

少し判定を厳しくすると

> (define (good-enough? guess x)
    (< (abs (- (square guess) x)) 0.00000001))
> (mys 2)
1.4142135623746899

大分精度が上がった。

そんじゃ、収束速度を調べてみる。sqrt-iterの冒頭に (print guess "n")を突っ込んで、何回呼ばれたかモニター。

> (mys 2)
1.
1.5
1.4166666666666665
1.4142156862745097
1.4142135623746899
1.4142135623746899
> (mys 2000000)
1.
1000000.5
500001.2499995
250002.62499475002
125005.31245537553
62510.65588770507
31271.325216410296
15667.640786075117
7897.646211720405
4075.4431105036065
2283.093646797556
1579.5490058334221
1422.8666047172292
1414.2398737439435
1414.2135626178515
1414.213562373095
1414.213562373095
> (mys 0.0002)
1.
.5001
.2502499600079984
.125524580467459
.0635589469493114
.03335281609280434
.019674655622311493
.014920008896897232
.014162413320389438
.014142150140500532
.014142150140500532