入れ子のかっこ

Gosperのハック

定数ステップで次を求めるかっこの足跡のアルゴリズムは, Gosperのハックにヒントがあるらしい. GosperのハックはMITのAIラボのHAKMEMのitem 175にあり, To get the next higher number with the same number of 1 bits のプログラムである. 元はPDP10の機械語だが, TAOCPのex7.1.3--20に通常の書き方がしてある. それをまたSchemeにすると

(define (gosper x)
 (let* ((u (iand x (- x)))
        (v (+ x u)))
  (+  v (irsh (quotient (ixor v x) u) 2))))

iand, ixor, irshは符号なし整数値を二進数とみてビット演算する. 途中結果も含めて計算の進行状況を見るとこうなる.

   x      u      v   xor x    / u     y
000111 000001 001000 001111 001111 001011
001011 000001 001100 000111 000111 001101
001101 000001 001110 000011 000011 001110
001110 000010 010000 011110 001111 010011
010011 000001 010100 000111 000111 010101
010101 000001 010110 000011 000011 010110
010110 000010 011000 001110 000111 011001
011001 000001 011010 000011 000011 011010
011010 000010 011100 000110 000011 011100
011100 000100 100000 111100 001111 100011
100011 000001 100100 000111 000111 100101
100101 000001 100110 000011 000011 100110
100110 000010 101000 001110 000111 101001
101001 000001 101010 000011 000011 101010
101010 000010 101100 000110 000011 101100
101100 000100 110000 011100 000111 110001

これをみて理由を考える. まずxがα01^a0^bであったとする. ただしa≥1, b≥0である. (1^aは1がa個並んでいること.) 上からの4行についていえば
x=000111, α=00, a=3, b=0
x=001011, α=001, a=2, b=0
x=001101, α=0011, a=1, b=0
x=001110, α=0, a=3, b=1
である. uは右端の1だから u=10^b. vはxの右端の1に1を足すから, a桁連続している1が0になり, 繰上げが出る. v=α10^a+b. xとvをxorするとαがキャンセルされて, 1^a+10^bになる. uで割ると0^bがなくなり, 2ビット右シフトすると, 1^a+1が1^a-1になり, 結局 α10^b+11^a-1が得られる.
xのαの右にあった0は1に変り, bが1増え, aが1減り, 場所が交代した.
結果的に1の数は変らない.

Gosparハックの逆関数もある.(TAOCPex7.1.3--21)

(define (igosper y)
 (let* ((t (+ y 1)) (u (ixor t y)) (v (iand t y)))
  (- v (quotient (iand v (- v)) (+ u 1)))))

次にように計算が進行する.

   y      t      u      v   and -v   u+1    x
111000 111001 000001 111000 001000 000010 110100
110100 110101 000001 110100 000100 000010 110010
110010 110011 000001 110010 000010 000010 110001
110001 110010 000011 110000 010000 000100 101100
101100 101101 000001 101100 000100 000010 101010
101010 101011 000001 101010 000010 000010 101001
101001 101010 000011 101000 001000 000100 100110
100110 100111 000001 100110 000010 000010 100101
100101 100110 000011 100100 000100 000100 100011
100011 100100 000111 100000 100000 001000 011100
011100 011101 000001 011100 000100 000010 011010
011010 011011 000001 011010 000010 000010 011001
011001 011010 000011 011000 001000 000100 010110
010110 010111 000001 010110 000010 000010 010101
010101 010110 000011 010100 000100 000100 010011
010011 010100 000111 010000 010000 001000 001110

方法としては, yの右端が0の時, 最も右の`10'を'01'にする. yの右端が1の塊の時, その左の0を挟んだ左の1から右を01と1の塊にし, その右に0を詰める.
これをifで分岐せずに実行するのが味噌であるが, 説明は場合を分ける.

yの右端が0の時, y=α10^b(b≥1), t=α10^b-11, u=1, v=y, v and -v=10^b, quotient (v and -v) (+ u 1)=10^b-1; x=α010^b-1.

yの右端が1の時, 1が全部右に寄ったら終りなので, 右の1の塊の左に0があり, その左に1があるとする. つまりy=α10^a1^b(a,b≥1)とする. t=α10^a-110^b, u=1^b+1, v=α10^a+b, v&-v=10^a+b, u+1=10^b+1 quotient (v&-v)(u+1)=10^a-1, これをvから引くとα10^a+b-10^a-1=α01^b+10^a-1.

入れ子のかっこ

Lisp屋には違和感が全くない, 入れ子のかっこの話題だ. TAOCPにかっこの足跡(parenthesis trace)という話がある. (-sesと複数ではなく, -sisと単数なのはなぜか.)
例えば, 左かっこ4個, 右かっこ4個で, 正当なかっこの組み合せを作ると, ()()()(), ()()(()), ..., (((())))の14通りが出来る. n=4で14通りになるというのは, 2009年8月のブログ, 「投票数」に書いた通り.
₈C₄-₈C₃=70-56=14だ.

さて, かっこの足跡は左かっこを0, 右かっこを1で表した二進数である. 従って, ()()()()は01010101, (((())))は00001111となる. この方法で辞書式順で書くと

色の線で囲ったのは, 同じパターンが見えるところだ.

かっこ構造を保ったまま, 辞書式順で次の構造に移るにはどうするか.　二進法の表記では, 右から0と1を数えながら左へ`01'を探す. (0と1の数を#0, #1と書こう.)　ただしこれまで通過した#0<#1でなければならない. そういう`01'を見つけたらそれを`10'と交換し, その右に#0だけ0を並べ, さらにその右に#1だけ1を並べる. これは辞書式順で最小の数を作るためである.

テストするには, 左から探したいので, 左右逆転のリストを使い,

(define x '(1 1 1 1 0 0 0 0))

(define (next x)
 (let ((z 0) (o 0))
  (define (dup n a) (if (= n 0) '() 
   (cons a (dup (- n 1) a))))
  (define (nx x)
   (if x (if (and (= (car x) 1) (= (cadr x) 0) (< z o))
           (append (dup o 1) (dup z 0) '(0 1) (cddr x))
           (begin (if (= (car x) 0) 
              (set! z (+ z 1)) (set! o (+ o 1)))
            (nx (cdr x)))) 'ok))
  (nx x)))

で(set! x (next x))を, x=okになるまで次々と実行する.

実はこれを一気に, 定数ステップで計算するアルゴリズムがあるらしい. これこそプログラムハックである.



μ₀は...010101というパターンの二進数である.

その計算の進行の様子を下に示す.

上の段の左端は説明用の行数である. 左上のxで, 赤い`01'は交換する場所である. 次のozは, それより右の1と0の数. xのビット位置を右から0,1,2,...と数える. #0=#1となる可能性があるのは, 奇数番目の位置にある0なので, 3行目の右端の`01'の0, 8行目の`0101'の0が要注意だ. 奇数番目を取り出すべく, μ₀を利用する.

tはxの右から見て最初の奇数桁目の1の右を0にする. uは反対に最初の奇数桁目の1の左を0に, 右を1にする. つまり右の1は要注意の0を隠す. この1の数は偶数だが, (#0+1)*2であることにも注意. vはuとxの∨でxの交換すべき`01'の右がすべて1になった. 交換すべき`01'の左は現状のまま. これに1を加えれば, `01'が`10'になり, 右は#0+#1個の0になる. 次のv∧¬wは右から#0+#1+1個の1なので, #0+1桁右シフトすれば, 右端に#1個の1が得られる. それが√(u+1)で割る理由である.

すばらしい!

細線化アルゴリズム

細線化アルゴリズムがどう働くか, 「和」の字を細くしただけではよく分からない. そこでまた実験をした. 横40ピクセル, 縦20ピクセルの黒い長方形を作り, このアルゴリズムで細線化した. その様子を下の図に示す.

玉ねぎの皮を剥くように「骨から肉が削がれる」様子が見える. 第1回目は1と書いた矢印の先の赤のピクセルたちが消える. 上から下へ各行を処理, 行内では左から右へ処理するとしよう.
f(x_NW,x_N,x_NE,x_W,x,x_E,x_SW,x_S,x_SE)=x∧¬g(x_NW,...,x_W,x_S,...,x_SE)
もともと0のピクセルは, fの式に x∧ があるから0のままである. この長方形の左上の隅で1のピクセルに最初に出会う. この時使うgのパターンは下の図に左端のものだ.


gの前に ¬ があるから, このパターンの時は要らないと読む. 従って左上隅は消える. 次に上の縁を処理するが, これはgのパターンの左から2番目で消える. パターンの下の数字は, 前回のブログの図のどこにあったかを, 左端を0として示す. このようにして, 上の縁, 右上隅, 右の縁, 右下隅の赤い地帯が消えていく.

2回目はgのパターンを180度廻転するから, 外周の青で示す左と下が無くなる. 3回目はオレンジ色の部分が消え, 4回目に緑が消える. こうして20行あった横の列は, 19回で骨だけになるのであった.

楕円でも事情は同様であった.


この楕円のデータは, TAOCPに楕円曲線の塗り潰しの例にあったものを使った. TAOCPによると, 楕円や円など, 楕円曲線の塗り潰しはことの他簡単という. そのうちやってみよう.

細線化アルゴリズム

TAOCP 7.1.3項の式159のGuo & Hallの細線化アルゴリズムがある. Life Gameのように, ピクセルのキング近隣(セルオートマトンでは,
Moore Neighbourhoodという. 自分と周りの8隣り)を見, 次の時刻の状態を
f(x_NW,x_N,x_NE,x_W,x,x_E,x_SW,x_S,x_SE)=x∧¬g(x_NW,...,x_W,x_S,...,x_SE)
で決める. ただし, gは, 近隣が

の時, 1とする.
上の操作を奇数回目とすると, 次の偶数回目は, gのパターンの向きを180度廻転して使う. それで2回連続で変化がなければ停止する.

本当かなぁ, と思い, 実装, テストした. 和田研フォントの「和」の字を200×200のグリッドにしてやってみた. その結果が次である.



まずまずはうまく行った. 演習問題に, M行N列の黒四角を細線化するとどうなるかというのがあり, それもやってみた. 30行40列だと, 中央あたりに11個の横1列が出来て終わりであった. 当たり前か.

両替問題

両替問題の続きである. 前回の話題はSICPの再帰関数の解法であった. コンピュータの数学やPolyaの解法は, 母関数による. コンピュータの数学は, 「この問題はPolyaが母関数を用いて解くよい教材であることを示して, 一躍有名になった」と書く.

母関数

数列を閉じた式で表す摩可不思議な方法である. Fibonacci数は0,1,1,2,3,5,8,...であるが, The Art of Computer Programmingの1.2.8の式(11)には, その母関数がG(x)=x/(1-x-x²)であると書いてある. (元はxでなくzだが) これをばか正直に割ってみると, 商のx^nの係数がfib(n)になる. つまり無限数列0,1,1,2,3,5,8,...を「閉じた式」にしたものが母関数である.



日本語は母関数だが, 英語ではgenerating functionという. なにかを産み出すという意味で, 母関数という訳語は悪くない. 統計には母集団があり, 円柱や円錐には母線があり, 分数には分母があり, 音声には母音がある. 母国語もある. 航空母艦や潜水母艦や捕鯨母船もあるが, 残念にも父の出番はあまりない.

「閉じた式」にするとなにがいいか. 例えば1,1,1,1,...に1,1,1,1,...を掛けると1,2,3,4,...になりそうである. 一方 1,1,1,1,...は閉じた形では1+x+x²+x³+...=1/(1-x)である. 従って1,1,1,1,...掛ける1,1,1,1,...は
1/(1-x)*1/(1-x)=1/(1-2x+x²)
これをばか正直に割ってみると
1+2x+3x²+4x³+...
が得られる.

私も真面目に読んだのではないが, TAOCPの1.2.9項や, コンピュータの数学の7章(ともに母関数)のところは, そのうちきちんと読みたい.

そのもっとも入門的な使い方が, Polyaの講義ノートにある両替問題なので, 簡単に説明しよう.

1セント硬貨で1セント, 2セント, ...にする方法はどれも1通りなので,
1+x+x²+x³+x⁴...
のように書く. 最初の1は1セント硬貨をまったく使わない場合である.

5セントでは
1+x⁵+x¹⁰+x¹⁵+x²⁰...
xやx²の項は, 作れないから係数が0で, 省略してある.

これを10セント, 25セント, 50セントで書き, (1)式のように掛け合わせれば, x¹⁰⁰の係数のところに, 両替法の数が出るわけである.

(2)は1セントの可能性の式を閉じた形にしたもので, (1)の因子をすべて閉じた式にしたのが(3)である. 最後は(8)が欲しいのだが, それを求めるために, 1セントだけの式(4), 1セントと5セントだけの式(5), などを作る. (4)の係数Aはすべて1だから, それと(5)とからBが(9), (10)のように得られ, 順にC, D, Eが(11), (12), (13)のように得られる.

早速Schemeのプログラムを書く.

(define (a n) (if (< n 0) 0 1))
(define (b n) (if (< n 0) 0 (+ (a n) (b (- n 5)))))
(define (c n) (if (< n 0) 0 (+ (b n) (c (- n 10)))))
(define (d n) (if (< n 0) 0 (+ (c n) (d (- n 25)))))
(define (e n) (if (< n 0) 0 (+ (d n) (e (- n 50)))))

(e 100) => 292

(a b c d e)の呼ばれた回数は(292 332 49 12 4)であった.

メモ化すると(21 40 24 8 4)に減った.
メモ化にはalistを使った.

(define (e n) (set! ec (+ ec 1))
 (let ((x (assoc n ea)))
  (if x (cadr x)
   (let ((r (if (< n 0) 0 (+ (d n) (e (- n 50))))))
    (set! ea (cons (list n r) ea)) r))))

のように書いてある.