入れ子のかっこ

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))))

のように書いてある.

両替問題

研究所の山本君がSICP(Structure and Interpretation of Computer Programs)の両替問題がなんとかといっていたので, 両替問題をおさらいしよう. 英語ではchangeというので, Graham, Knuth, Patashnik著, 有澤他訳の「コンピュータの数学(Concrete Mathematics)」では, 「釣銭の問題」(299ページ)となっている. 正しくは, 空港のCHANGEの窓口と同じで, 外国のお金から, または高額のお金からの両替である. 数学の問題としては, 高額のお金の崩し方の話題である.

この問題はいろいろなところに書いてある. 前述のSICPや, コンピュータの数学の他, G. PolyaのHow to Solve It(如何にして問題を解くか)にもあり, 私がたまたま持っている, PolyaとTarjanの1978年のStanfordでの講義ノート, D.R.Woods, Notes on Introductory Combinatoricd(STAN-CS-79-732)にも出ている.

SICPの記述はこうだ

50セント, 25セント, 10セント, 5 セント, 1セントがあるとして1ドルの両替にはどのくらいの場合があるか. つまり, 金額に対して両替の場合の数を計算する手続きは書けるだろうか.

この問題は再帰的手続きとして単純な解がある. 使える硬貨をある順に並べたとしよう. すると次の関係が成り立つ:

n種類の硬貨を使う, 金額aの両替の場合の数は:
• (a)最初の種類の硬貨以外を使う, 金額aの両替の場合の数, 足す
• (b)dを最初の硬貨の額面金額[denomination]として, n種類の硬貨を使う, 金額a-dの両替の場合の数

つまり1セントと5セントの2種類の硬貨で10セントにするには, 最初の種類の硬貨は5セントであり, 5セント以外の硬貨, つまり1セントで10セントにする場合の数(下の図で赤枠の場合)と, 5セントと1セントで10セントから5セント引いた5セントにする場合の数(青枠の場合)の和である.

また再帰の終わり

• (c)aがちょうど0なら, 両替の場合の数は1
• (d)aが0より少なければ, 両替の場合の数は0
• (e)nが0なら, 両替の場合の数は0
も考えなければならず, 結局プログラムは以下のようになる.

(define (count-change amount)
  (cc amount 5))

(define (cc amount kinds-of-coins)
  (cond ((= amount 0) 1)                ;(c)
        ((or (< amount 0) (= kinds-of-coins 0)) 0);(d)と(e)
        (else (+ (cc amount             ;(a)の場合
                     (- kinds-of-coins 1))
                 (cc (- amount          ;(b)の場合
                        (first-denomination kinds-of-coins))
                     kinds-of-coins)))))

(define (first-denomination kinds-of-coins)
  (cond ((= kinds-of-coins 1) 1)
        ((= kinds-of-coins 2) 5)
        ((= kinds-of-coins 3) 10)
        ((= kinds-of-coins 4) 25)
        ((= kinds-of-coins 5) 50)))

10セントの両替で, ccの呼ばれ方を図にすると下のようになる. 頂点の(10 2)はamountが10, kinds-of-coinsが2でccが呼ばれるの意. 左下(10 1)への枝は(a)の下請け, 右下(5 2)への枝は(b)の下請けであり, さらに下請けが次々と呼ばれることを示す. さらに右下方向へ延びる枝が, 図には全体で3本見えるが, この最後が(c)の1通りを返してきて, 結局35回呼ばれ, 3通りになる.


このプログラムの, 100セントの両替では, 関数ccは 15499回呼ばれた. その解析がこの表である.


表の横はamount, 縦はkinds-of-coinsで, 横の見出しの「--」はその間の数, つまり0と5の間なら, 1,2,3,4のことである. amountは硬貨3種類のときに-5で3回呼ばれる. これは(cc 5 3)の時で, 5セントから10セントを引こうとして生じる.

この表の数値を全部足すと, 15499になる. また数値の入っている場所は250個所ある. つまり250通りで呼ばれるのである. 赤枠内のamount=0の欄を足すと, 解の292が得られる.

メモ化

これとかFibonacci数の計算のプログラムのような再帰プログラムでは, 同じ引数で関数を何回も呼ぶことによくなる. 上の計算でも引数の異なる呼出しは, 250回である. こういう時には, メモ化という手法が有効である.

小学生の頃, 長い割り算(商の桁数が多い)をするなら, まず除数の2倍, 3倍,..., 9倍の表を作ってから割り算の取り掛かったが, メモ化は表を新しい引数の組で関数を計算するたびに作るのである.

SICPでは, 3.3.3節の最後にメモ化の話題がある. それに倣って, memo-ccを作ってみた.

(define (make-table)  ;和訳SICP 159ページにあるmake-table
  (let ((local-table (list '*table*)))
    (define (lookup key-1 key-2)
      (let ((subtable (assoc key-1 (cdr local-table))))
        (if subtable
            (let ((record (assoc key-2 (cdr subtable))))
              (if record
                  (cdr record)
                  false))
            false)))
    (define (insert! key-1 key-2 value)
      (let ((subtable (assoc key-1 (cdr local-table))))
        (if subtable
            (let ((record (assoc key-2 (cdr subtable))))
              (if record
                  (set-cdr! record value)
                  (set-cdr! subtable
                            (cons (cons key-2 value)
                                  (cdr subtable)))))
            (set-cdr! local-table
                      (cons (list key-1
                                  (cons key-2 value))
                            (cdr local-table)))))
      'ok)
    (define (dispatch m)
      (cond ((eq? m 'lookup-proc) lookup)
            ((eq? m 'insert-proc!) insert!)
            (else (error "Unknown operation -- TABLE" m))))
    dispatch))

(define (memoize f) ;和訳SICP 160ページのmemoizeを2次元化
  (let ((table (make-table)))
    (define get (table 'lookup-proc))
    (define put (table 'insert-proc!))
    (lambda (x y)
      (let ((previously-computed-result (get x y)))
        (or previously-computed-result
            (let ((result (f x y)))
              (put x y result)
              result))))))


(define memo-cc       ;和訳SICP 22ページのccをメモ化
  (memoize (lambda (a k)
    (set! cccount (+ cccount 1))
    (cond ((= a 0) 1)
          ((or (< a 0) (= k 0)) 0)
          (else (+ (memo-cc a (- k 1))
              (memo-cc (- a (first-denomination k)) k)))))))

これでやってみると, memo-ccが呼ばれたのはちょうど250回であった.

微分解析機

遊星歯車の動き方のアニメーションを作りたいと思っていたが, やっとなんとかなったので, 公開することにした. 試作の順に番号をつけたので, planetary2とplanetary4とがある.

planetary2では


またplanetary4では


のような絵が得られる.

とりあえず動きを見るために作ったのがplanetary2で, 歯車を円で近似した. しかし, 歯車の廻るところも見たくなり, planetary4を書いた. ただ歯車は早く廻転すると, ストロボ効果で, 逆転のように見えるので, 結構遅く廻転させている.

下にある赤, 緑, 橙, 青の四角をクリックすると, いずれかの行動をする. そのパターンが

である. 赤線は太陽歯車の廻転角, 緑線は内歯車の廻転角, 青線は衛星キャリヤの廻転角である.