newLISPを改造する

とある方から、

 話題にしていた文献について添付しました。
 
 この中でfigure 11 (a),(b)は、回転していないとの表現がありますが、
 理解に苦しんでいます。

こんな質問を寄せられた。問題の文献は、添付だったけど、探したら、 http://www.ieee.li/pdf/essay/quadrature_signals.pdf にあった。有名な先生のエッセィだそうだ。

ふむ、z軸が時間なら、時と供にベクトルが回転してくってのは認められるけど、この図の z軸は、周波数。周波数は回転しない。と禅問答。

まてよ、各周波数毎に、その成分が複素面に投影されるって事を、この図は言ってるんだな。 と言う事は、時間領域から周波数領域へのワープ。フーリエ変換の出番だな。

丁度、今newLISPをやってて、その中にさりげなく、fft,ifftが仕込まれているよ。

> (fft '(1 2 3 4))
((10.000000 0.000000) (-2.000000 -2.000000) (-2.000000 0.000000) (-2.000000 2.000000))

等差数列のデータを周波数領域へワープさせる。結果は、4つの複素数で返ってきた。 一番最初のデータは直流成分、次は、周期4の周波数成分。オイラー平面だと、第3象限に 現れるのか。

これをパワースペクトルに変換してみる。

>
(define (ps v)
  (sqrt (add (pow (v 0) 2) (pow (v 1) 2))))

(lambda (v) (sqrt (add (pow (v 0) 2) (pow (v 1) 2))))
> (map (fn (z) (ps z)) (fft '(1 2 3 4)))
(10.000000 2.828427 2.000000 2.828427)

実際の値を求めるなら、データの個数で(この場合は4)で割ればいいんだな。

そんじゃ、FFTの実行スピードは? 簡単にベンチしてみた。

> (time (fft (sequence 1 256)))
0.206000
> (time (fft (sequence 1 1024)))
1.236000
> (time (fft (sequence 1 4096)))
3.483000

かなり良い線行ってるな。ノイズが乗った信号をFFTして、周波数ドメインへ移し、そこで フィルターをかけてから時間ドメインへ戻すって事も出来そう。 こんなやりかたも別目的で使われる事がある。音データとか画像データの圧縮。mp3だとか jpegなんてのは、この技術を使ってる。

折角なんで、どんなソースなってるか探してみると、核になる部分は、nl_math.cの中で fastFourierTransformって 関数になってた。

stepperを使い易く改造する

newLISPをキーにいろいろしてたら、 様々な LISP なんてのに引っかかった。みんな大好きオレオレLISP。そして昔流行った、基礎文法最速 マスター編にも、 newLISP基礎文法最速マスター が、登録してあったぞ。この方も、オレオレのお仲間なんでしょうな。 オイラーも、遅ればせながら、お仲間に入れて貰おう。

最初の話題は、前回やったnewLISPのdebugger。debuggerと言いつつ、実態はstepperだけどね。 stepさせるかnextさせるか、その都度入力するの、とーーーーても面倒。RETURNを叩いたら、 前回に同じってなれば、楽できるな。世の中のgdbもそうなってるし。

改造しましょ。ソースが有るなら改造を楽しむのが八苦の始まりですよ。

まずはどのファイルを改造するか当たりを付ける。何、nl-debug.cって名は体を表しているんで、 悩む事はありません。次はプロンプトを目当てに絞ります。そうすると、fgetsを 使って、コマンド入力を受け取ってる部分が付近にありました。

ロジックを読みます。改行が有ったら、その前の文字を見て、s,n,c,qかで分岐してました。 ならば、改行のみの時に、s,n,c,qのどれかを内部的に入力してあげるのが簡単そう。 前のコマンドを覚えておけば、そういう事が出来るね。

そうやって書いたのが下記。ついでに、S式の範囲を明確にするエスケープ文字に変更しておきました。

sakae@uB:~/newlisp-10.6.2$ diff -u nl-debug.c.orig nl-debug.c
--- nl-debug.c.orig     2015-07-22 05:05:11.317825905 +0900
+++ nl-debug.c  2015-07-22 06:02:56.592589576 +0900
@@ -43,11 +43,12 @@
 int debugStackIdx = 0;
 UINT *debugStack;

-char debugPreStr[8] = "#";
-char debugPostStr[8] = "#";
+char debugPreStr[8] = "\e[7m";
+char debugPostStr[8] = "\e[0m";
 CELL *debugPrintCell = NULL;
+char prvcmd;

-char headerStr[16] = "\n-----\n\n";
+char headerStr[16] = "\n-----\n";
 char footerStr[32] = " s|tep n|ext c|ont q|uit > ";

 #define pushDebugStack(A) (*(debugStack + debugStackIdx++) = (UINT)(A))
@@ -300,16 +301,22 @@
        fatalError (ERR_IO_ERROR, 0, 0);

       /* client and server could have different line-termination */
+      if (*(command) == '\n' || *(command) == '\r'){
+         *command = prvcmd;
+         *(command + 1) = '\n';
+      }
       if (*(command + 1) == '\n' || *(command + 1) == '\r')
        {
          if (*command == 'n')
            {
+             prvcmd = 'n';
              traceFlag |= TRACE_DEBUG_NEXT;
              currentLevel = recursionCount;
              break;
            }
          if (*command == 's')
            {
+             prvcmd = 's';
              traceFlag &= ~TRACE_DEBUG_NEXT;
              break;
            }

とまあ、簡単に出来たように書いてますが、裏には八苦がありました。この際ですから、 恥ずかしい事を書いておきます。前のコマンドを覚えるタイミングを、4種のコマンド 解析が終わった後でやっていた。望み通りの動きにならない。

ならばemacs上からgdbを使って、動いているnewlispにattachしようと したけど、emacsがHungしちゃって駄目。はてなと思ってgdbのバージョンを調べたら、 システム備え付けで6.3が入っていた。今は、gdb7.9.1って時代なのに。更新してないのは、 カーネルdebugをする為に特製のpatchが当たっているから? FreeBSDも時代の遺物みたいな やつが使われていたぞ。

で、新しいのを入れたらどうよ? やってみた。/usr/local/binの下にpkg物は入るんだけど、 gdbなんてコマンドは無い。普通はgdb-7.9.1みたいな名前で入るはずだけどね。 pkg_info -L gdb で内容物を確認したら、egdb ですって、これでようやくemacsから M-x gdbして

Run gdb (like this): egdb -i=mi newlisp 12076

で動くようになった。(12076は、別端末で動かしてるnewlispのpid)

動作確認は、例によってfact。それはいいんだけど、factの引数名に n なんてのを 使っちゃったもので、stepperのコマンドとかぶってしまって変数の値確認が出来ない。 s,n,c,qは、変数名にしないようにしましょう!

大きな数とか

折角factを書いたので、大きな数の計算も出来るかやってみる。

> (fact 50)
-3258495067890909184

これはBugか? マニュアル見たら、大きな数は、サフィックスにLを付けるとな。

> (fact 50L)
30414093201713378043612608166064768844377641568960512000000000000L

数の扱いで16進数とかも受けてくれる。

> 0xff
255
> 010
8
> 0b1000000000000
4096

それじゃ、逆演算は? 英語でもそうだけど、Highって単語を覚えたら、その逆も一緒に 覚えると効率的ですからね。取りあえず、16進数が得られればいいな。

> (format "0x%x" 255)
"0xff"
> (format "0x%X" 4095)
"0xFFF"

dump

マニュアルをつらつら見てたら、CELLの内容を分かり易く表示してくれる関数が載っていた。

以前、Lispのご本尊様は? ってんでソースを紐解いた時に出てきてたやつ。違和感が 有ったんだよな。普通のLispとかだと、car,cdr相当が用意されてるんだけど、そんな項は 無かった。マニュアルによると、5つのデータが得られるとか。意味は、

offset description(詳細)
0      newLIPセルのメモリ・アドレス
1       cell->type:
          major/minor 型、詳細は newlisp.h を見てください。
2       cell->next:
          リスト・ポインタへのリンク
3       cell->aux:
           文字列の長さ+1か、
           64ビット整数のロー(リトル・エンディアン)またはハイ(ビッグ・エンデ
           ィアン)ワードか、IEEE 754 倍精度浮動書数点数のロー・ワード
4       cell->contents:
           文字列/シンボルのアドレスか、
           64ビット整数のハイ(リトル・エンディアン)またはロー(ビッグ・エンデ
          ィアン)ワードか、 IEEE 754 倍精度浮動書数点数のハイ・ワード

らしい。

> (dump first)
(2087043904 263 2087043072 976733431 439912496)
> (dump 123)
(2087049744 898 2087043072 123 0)

これ、10進数ぽいな。dumpに10進数とは、素人臭いな。ええい、こうしてやる。

> (map (fn (v) (format "%lX" v)) (dump 9223372036854775807))
("7D440A00" "382" "7D43F000" "FFFFFFFF" "7FFFFFFF")
> (map (fn (v) (format "%lX" v)) (dump -9223372036854775808))
("7D4409F0" "382" "7D43F000" "0" "80000000")
> (map (fn (v) (format "%lX" v)) (dump 9223372036854775807000L))
("7D440A70" "582" "7D43F000" "4" "7D7A1D20")

マニュアルから大きな数を拾ってみた。int64なのね。

マニュアルには、関連項目もあって、cpymem address pack unpack あたりを見ておくと 幸せになれるらしい。(それとも八苦か?)

ならば、今後の為に、ソースをオイラー流に整えておきましょ。

[ob: newlisp-10.6.2]$ for f in *.[ch]
> do
> indent -bsd -i4 $f
> done

bsdスタイルで、4文字タブに整形しました。で、dump関数に引数を与えないで実行すると、 システムに収監されてるオブジェクトをくまなく表示してくれるそうです。

> (dump)
address=7C65C000 type=256 contents=nil
address=7C65C010 type=257 contents=true
address=7C65C020 type=6 contents=MAIN
address=7C65C040 type=263 contents=while@1A387550
address=7C65C050 type=263 contents=until@1A387510
   :
address=7C65DD60 type=69 contents=-
address=7C65DD70 type=69 contents=n
address=7C65DD80 type=898 contents=1
address=7C65DD90 type=316 contents=(lambda (n)
 (if (zero? n)
  1
  (* n (fact (- n 1)))))
address=7C65DE40 type=69 contents=args
address=7C65DE50 type=316 contents=(lambda () (cons (context) (args)))
true

何やら、先ほど書いたfactの一部が見えますなあ。それはそれで良いけど、どうもtypeが 10進数っぽい。これは頂けませんね。ソース見ましょ。

newlisp.cの中のp_dump関数にやってる事が書いてありました。最初の部分は、引数を与えた場合の処理。 次の部分はcellMemoryからobjを辿って表示してきます。その中で、typeも見てるんだけど、 10表示指定になってたので、修正。

canvas.lspを読み込んで、どんなタイプのセルが出来るか確認してみたい。newlispのログ 指定、-Lは、こういう場合無力なのね。下記を実行して、画面のdumpを取ってみた。

[ob: pc]$ newlisp -e '(module "canvas.lsp") (dump) (exit)' >LOG

どのぐらいのCELLが使われるかを調べておく。

> (sys-info)
(439 268435456 411 1 0 2048 0 296 10602 130)
> (module "canvas.lsp")
MAIN
> (sys-info)
(1424 268435456 514 1 0 2048 0 296 10602 130)
> (- 1424 439)
985

ふむ、985個のCELLを消費したとな。そんじゃ、取得したログを加工してtypeの部分だけを抜き出してみる。 左の数値は、使われている個数です。

[ob: pc]$ cut -d' ' -f2 LOG | fgrep type= | sort | uniq -c | sort -nr
 488 type=45     ; CELL_SYMBOL
 378 type=107    ; CELL_PRIMITIVE
 299 type=3b     ; CELL_EXPRESSION
  99 type=104    ; CELL_STRING
  74 type=100    ; CELL_NIL
  55 type=382    ; CELL_INT (int or long or int64 or bigint)  
  43 type=13c    ; CELL_LAMBDA
  26 type=1a     ; CELL_QUOTE
   7 type=6      ; CELL_CONTEXT
   4 type=183    ; CELL_FLOAT
   1 type=182    ; CELL_INT
   1 type=101    ; CELL_TRUE

newlisp.hのcell typesから、オイラーの目で拾って(下4Bitだけ見て)、転記してみました。 上位のビットは多種のフラグになってるんで、intが2箇所に出てるけど、詳細は マシンになって解読して下さい。

そんじゃ、先に使った

(define (fact n)
  (if (zero? n)
      1
      (* n (fact (- n 1)))))

をロードした時、どんな風にセルが使われているか調べてみるか。

[ob: pc]$ newlisp -e '(dump) (exit)' > BASE
[ob: pc]$ newlisp -e '(load "fact.lsp") (dump) (exit)' > ADD

こんな風に、初期状態のCELLの様子とfact.lspをロードした後のCELLの様子を採取。 後は、このログファイルのdiffを取れば、ロードした時の差分が見られるはず。

ありゃりゃ、全部違うって言ってきたぞ。これはいかなる事が起こってる?

CELLの割り当ての源流を探る

普通のlispだと、起動時にがーとメモリーを取ってきて、未使用CELLの列を作る。 CELLが必要になった時、未使用CELLの列から一つ貰ってきて使うって事をする。 newLISPも同じ事をやってるのだろう。

CELLを消費する関数(例えばcons)を見ると、getCellがそれっぽい。んで、起動時から どんな風に呼ばれるか、gdbで追ってみる。

(gdb) bt
#0  _thread_sys_mmap (addr=0x0, len=65536, prot=3, flags=4098, fd=-1, offset=1969704960) at /usr/src/lib/libc/sys/mmap.c:53
#1  0x011ad70d in map (d=<optimized out>, sz=65536, zero_fill=0) at /usr/src/lib/libc/stdlib/malloc.c:458
#2  0x011ae1ac in omalloc (sz=65536, zero_fill=0, f=0x0) at /usr/src/lib/libc/stdlib/malloc.c:1074
#3  0x011af35a in malloc (size=65536) at /usr/src/lib/libc/stdlib/malloc.c:1173
#4  0x151d9462 in allocMemory (nbytes=65536) at newlisp.c:2355
#5  0x151d9573 in allocBlock () at newlisp.c:2268
#6  0x151d9b07 in getCell (type=256) at newlisp.c:2083
#7  0x151e3a2b in initialize () at newlisp.c:1341
#8  0x151e4ecd in main (argc=1, argv=0xcfbfcfe4) at newlisp.c:726

newlisp内のallocMemoryからmallocのサービスを要求してんだけど、mallocの中まで下りて 行けたぞ。普通のunixだと、libの中までは入って行けないのに、OpenBSDはサービス精神 旺盛だなあ。

mallocの中では、普通、sbrkとかでヒープを拡張して、拡張された部分を割り当てるってのが、C語の バイブルに書いてあるんだけど、どうやら違った。mmapってのを利用してるっぽい。

そこで、mallocはどう実現されてるか?なんてのを、 メモリー管理の内側 を参照してみたよ。ガベコレまで紹介されてて、一時あの人の本が売れてたな。 ああ、ウィキペデアのmallocも参考になるな。

OpenBSDはセキュリティに気を使って、毎回使うアドレスを変えてるんだな。

#define MMAP(sz)        mmap(NULL, (size_t)(sz), PROT_READ | PROT_WRITE, \
    MAP_ANON | MAP_PRIVATE, -1, (off_t) 0)

egdbでは、OSとの境界まで行けたけど、そこまで行かなくても、肝は、malloc.cに定義 されてたマクロだな。fdが -1 って、どんなファイル?普通に考えれば、無効なfdに相当 するな。ファイルにmapするのを諦めて、ただメモリーを確保するって事だろうね。 OS領域まで偵察を出して、 uvm/uvm_map.cを調べてみたけど、明確に成らなかった。しかしmanを注意深く読むと、fdが-1ってのを 特別扱いしてるよと暗に示されている。読みがまだまだ浅いな。

まあいい、newlispの動きとしては、CELLの要求が有った時に4096個分のセルに相当する メモリーを貰ってきて、それをCELL用に地均しして使うとな。

他のOSでの挙動はどうか? 調べてみたらウブもfedoraも起動の度にmallocは違うアドレスを 返してきた。まあ、ディストロが異なっても、glibcは一緒だから同じ挙動か。

FreeBSDはどうか? mallocは常に同一のアドレスを返してきた。これでやっとdiff出来る。

@@ -406,10 +406,10 @@
 address=2884F950 type=100 contents=nil
 address=2884F960 type=182 contents=0
 address=2884F970 type=3b contents=()
-address=2884F980 type=3b contents=("newlisp" "-e" "(dump) (exit)")
+address=2884F980 type=3b contents=("newlisp" "-e" "(load \"fact.lsp\") (dump) (exit)")
 address=2884F990 type=104 contents="newlisp"
 address=2884F9A0 type=104 contents="-e"
-address=2884F9B0 type=104 contents="(dump) (exit)"
+address=2884F9B0 type=104 contents="(load \"fact.lsp\") (dump) (exit)"
 address=2884F9C0 type=1a contents='(set
  (global 'module)
  (lambda ($x) (load (append (env "NEWLISPDIR") "/modules/" $x))))
@@ -492,5 +492,29 @@
 address=2884FE30 type=13c contents=(lambda () (cons (context) (args)))
 address=2884FE40 type=1a contents='(dump)
 address=2884FE50 type=3b contents=(dump)
-address=2884FE60 type=45 contents=dump
+address=2884FE70 type=45 contents=dump
+address=2884FEF0 type=100 contents=nil
+address=2884FFF0 type=45 contents=n
+address=28850000 type=3b contents=(n)
+address=28850010 type=3b contents=(if (zero? n)
+ 1
+ (* n (fact (- n 1))))
+address=28850020 type=45 contents=if
+address=28850030 type=3b contents=(zero? n)
+address=28850040 type=45 contents=zero?
+address=28850050 type=45 contents=n
+address=28850060 type=382 contents=1
+address=28850070 type=3b contents=(* n (fact (- n 1)))
+address=28850080 type=45 contents=*
+address=28850090 type=45 contents=n
+address=288500A0 type=3b contents=(fact (- n 1))
+address=288500B0 type=45 contents=fact
+address=288500C0 type=3b contents=(- n 1)
+address=288500D0 type=45 contents=-
+address=288500E0 type=45 contents=n
+address=288500F0 type=382 contents=1
+address=28850100 type=13c contents=(lambda (n)
+ (if (zero? n)
+  1
+  (* n (fact (- n 1)))))
 true

セルの消費状況がよく分かって、興味深い。違う種類のOSを入れていると、時には 挙動の違いに悩まされる事があるけど、こういう場合は便利だね。