細線化アルゴリズム

細線化アルゴリズムがどう働くか, 「和」の字を細くしただけではよく分からない. そこでまた実験をした. 横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を書いた. ただ歯車は早く廻転すると, ストロボ効果で, 逆転のように見えるので, 結構遅く廻転させている.

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

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

2010

Knuth先生のホームページに
Did you know that 2009 is (1/2+3)*(4+5*(6*7+8*9))?
と書いてある. (年があけてからはHappy 2010 (i.e., MMIX++) to all!になった.)

今年は2010年. 今度はどうなるかと思ってプログラムを書いた.
結論からいうと

(- (- (+ 1 2) (* (* (* (- (- 3 4) 5) 6) 7) 8)) 9)
(* (/ (* 1 2) 3) (- (* (* (* (+ 4 5) 6) 7) 8) 9))
(+ (- 1 (* (+ 2 (* (* (- (- 3 4) 5) 6) 7)) 8)) 9)
(* (- (- 1 2) 3) (- (/ (+ 4 5) 6) (* (* 7 8) 9)))
(* (* (/ (* 1 2) 3) (+ (- 4 5) (* (* 6 7) 8))) 9)
(* (+ (+ 1 2) (* 3 4)) (+ (- 5 6) (* (+ 7 8) 9)))
(* (+ 1 (* 2 (+ 3 4))) (+ (- 5 6) (* (+ 7 8) 9)))
(* (+ (* 1 2) (* (- 3 (* 4 (- (/ 5 6) 7))) 8)) 9)
(* (* 1 (+ 2 (* (- 3 (* 4 (- (/ 5 6) 7))) 8))) 9)
(+ (+ 1 2) (* (- (* (+ (* 3 (+ 4 5)) 6) 7) 8) 9))
(+ (+ 1 2) (* (+ 3 (* 4 (+ (+ 5 (* 6 7)) 8))) 9))
(+ (+ 1 (* (* 2 (- (* (* 3 4) (+ 5 6)) 7)) 8)) 9)
(* (+ (+ 1 2) 3) (- (* 4 (+ (* 5 6) (* 7 8))) 9))
(* (* (* 1 2) 3) (- (* 4 (+ (* 5 6) (* 7 8))) 9))
(+ (+ 1 2) (* (+ 3 (* 4 (+ (- 5 6) (* 7 8)))) 9))
(* (+ (- (- 1 2) 3) (* 4 (+ (/ 5 6) (* 7 8)))) 9)
(+ (+ 1 2) (* 3 (+ (* (* 4 (+ 5 6)) (+ 7 8)) 9)))

の17通りが2010になる. 11番目と最後の加算と乗算だけの式は美しい.
規則としては, 数字は1から9までこの順に使う. 演算子は+, -, *, /だけである.

実はTAOCPの第4巻, 分冊4の7.2.1.6小項の演習問題122は,

100=1+2*3+4*5-6+7+8*9=(1+2-3-4)*(5-6-7-8-9)
                     =((1/((2+3)/4-5+6))*7+8)*9

などの例の後で, このような100の作り方はどれだけあるかというのである. 候補の式の数は, 10646016あり, 100になるのは1640通りあるそうだ.

100の美しい式は 1+2+3+4+5+6+7+8*9 である.

このような式で表せない最小の正の整数は, 2284だそうだから, まだ当分楽しめる.

上の17通りの全解探索に要した計算時間は, MacBookAirで50分弱であった.
MIT Schemeを使ったが, 分数が使えるので, こういう計算には便利である.


後記
折角プログラムを書いたので, 100になる式をやってみた. ところが1640個ではなく1641個が得られた. Knuth先生にメイルしたところ, 珍しく早々にメイルで返事が来た. もちろん先生の秘書から来る.
それによると
2004年に書いたプログラムにはミスがあり, 1641が正しい.
先頭に単項の-をつけると2010になる式には他に26通りある.
ということだ.

微分解析機

微分解析機の加算装置に, 車のデフ型でなく, 遊星歯車型の差動歯車を使う話を書いたが, 正月休みにもう少し定量的に理解したいと思った.

遊星歯車の動き方のアニメーションがここやここにある. でもこのくらいは自分で描きたい.

まず各歯車の動き方を調べた. 下の図で, 周囲の大きい円が外側の内歯車(歯が内側を向いている), その半分の半径の中央の円が太陽歯車Sである. 右の小型の円が衛星歯車で, 真横にあるのが最初の位置Pである.



太陽と衛星歯車は, AとBで, 衛星歯車と内歯車は, CとDで接している. 太陽歯車がA'まで廻転し(回転角をφ), 衛星キャリヤがP'まで来た時, キャリヤの廻転角をθとする. この時, 遊星歯車と内歯車はC"とD"で接している. その移動中に接した歯車は円弧DD"とC'C"でこの長さは等しい. 一方AとBはすでにA'とB'まで移動していて, 今はA"とB"で接している. この方の移動中に接した歯車は円弧A"A'とB"B'でこれも前の長さに等しい.
内歯車, 太陽歯車, 遊星歯車の半径をそれぞれ, R, r, ρ; B"P'B'の角をαとする.
DD"=Rθ=C'C"=ρα=B"B'=A"A'=r(φ-θ)
従ってキャリアの廻転角は θ=φr/(R+r);
衛星歯車の廻転角は α=θR/ρ=φRr/(R+r)ρとなる.

一方, 太陽歯車を固定し, 内歯車をψだけ廻転すると, キャリヤの廻転角θも同様に計算出来る.

θ=ψ*R/(R+r);
α=-ψ*Rr/(R+r)ρとなる.

この両者を合わせ, 太陽歯車の廻転角 φと内歯車の廻転角 ψから, キャリヤの廻転角θは
θ=φ*r/(R+r)+ψ*R/(R+r);
遊星歯車の廻転角αは
α=(φ-ψ)*Rr/ρ(R+r).

前回の歯車の図や, 上の図ではR=200,r=100,ρ=50なので,
θ=φ/3+2*ψ/3;
α=4*(φ-ψ)/3となる.

例えば, 太陽歯車を90度進めると, φ=90からθ=30度, α=120度(上の図). 次に内歯車を90度進めて追い付くと, ψ=90からθ=60度, α=-120度で, キャリヤも90度まで進み, 衛星歯車も最初の接していた両側の点で再び接するようになる.

通信機用のバリコンのつまみは, 太陽歯車に相当し, バリコンはキャリヤ軸についていたとすると, 1/3の減速になるのであった.