dancing links

KnuthのTAOCP, 7巻の分冊が出回り始めた. その最初の方にexact cover problemを解くのにdancing linksが適してると書いてある.

dancing linksについてはKnuthが楽しんで書いた論文があり, 講義のビデオも存在する. exact cover proble (日本語では「敷き詰め問題」か) は条件をならべた行列から解の候補を選び, それにより使えなくなった行や列を外し, またback trackするときは, 外した行や列をもとにもどす.

この行列はスパースなので, 上下左右に双方向リンクで実装すると, うまくいくというのがdancing linksの話の味噌である.

私としては面白そうなアルゴリズムを知ったら, さっそくコーディングして例題を走らせ, 体験によってアルゴリズムを理解, 記憶することが多い.

dancing linksについてもさっそくプログラムを書いた. googleで探すとjavaなんかで書いたプログラムも見つかるが, 自分で書くのが楽しいのでschemeで書くことにした.

Knuthの論文の最初の簡単な例題は, 直ぐにうまくいったが, 当然dancing linksを使えばうまく解けると書いてある他の例は, 結構てこずった.

8×8-2×2の60区画にペントミノのピースを置く問題である. ピースXを左上の(3,3)に置くというのに挑戦た. ぜんぜん解が出てこないのである. 当然プログラムがおかしいかといろいろ見てみたが, 合っているとしか思えない.

実は候補の列を選ぶのに, ブランチ数の少ないものを選ぶと書いてあったが, そうしなくても解が出るようにも読めたので, 残っている列の最左端から選んだのが失敗であった. あるとき思いついて, ブランチ数最小の列を選んでみたら, たちまち最初の解が出た. 超簡単な例では, 最左端でも問題なく解が得られたが, ペントミノ程度になると, ブランチ数を考慮する必要があった.

dancing linksの双方向リンクによる削除挿入は, 野下浩平君とその学生の一松宏君が8クイーン問題を解くのに考えたといわれている. しかし縦と横はすべてにクイーンを置くので, exact cover problemであるが, 斜め方向は半分くらいしか埋めないので, この辺の扱いは別になる. Knuthの例の論文をよく読むと, こういうのはgeneralized exact cover problemというので, 列にはprimaryとsecondaryを用意する. primaryの列名は双方向リンクで接続するが, secondaryの列名は自分自身で循環するリンクにすると書いてあった. そういう風にプログラムをなおすと, なるほどうまく行った.

数独にも適しているといわれるdancing linksである. 数独ではいくつかの解が与えられているので, 解を探し始める前に, 与えられた解について列や行を削除する(coverする)必要がある. 解に相当する行を探す手間が必要で, この辺はプログラムでかなり工夫が必要であった.

最後に, 上下左右に双方向リンクにしないでも, 実装できるのではないかとプログラムをしてみた. ノードはlispの対で, carに下へのリンク, cdrに右へのリンクを置き, 循環するようにリストを実装すると, たしかにこれでもうまく出来ることがわかった.

中山道

以前から少しずつ歩いていた中山道533Kmも遂に京都三条大橋へ辿り着いた. (詳しくいえば草津宿から大津宿を経て三条大橋までは旧東海道) 一緒に歩いてくれた多田君, 寺田君, 岩崎君に感謝したい. 実はまだ日本橋から本郷東大前を通り埼京線板橋駅まで8Kmが残っているが, これは車やバスで幾度も通り, 勝手知ったる道なので, もう済んだようなものでもある.

中山道は人気があり, 解説書やウェブページも多く, 歩こうと思えば問題は少ない. われわれは「ボクたちが歩く中山道」を手引きにした. また歩き出してから気づいたが, 道中処々に赤いシールが貼ってあり, 分岐点などで重要な道案内になっていた. 誰が貼ったのか知らないが有り難かった.

トラックがびゅんびゅん走る国道を歩くこともあったが, だいたいは静かな田舎道で, スローライフそのものである. それと引換えに山道では, 宿泊場所が左程ない; 食堂は国道を横切るときに探すしかない; という制約があることも判明した.

それにしても春秋のよい季節にのんびり歩けるのは最高であった. 昔の人は16日で歩き通したらしいから, 1日あたり33Km程度の速さである. これではのんびり歩く気分ではなかったを思われる. (中山道最大のイベント和宮東下でも25日) われわれは半日行程では10～15Km. 1日行程では20～25Kmであった.

昔と違うのは電車線と併走している場合, 草臥れると最寄駅から電車で宿泊地へ向かい, 翌日また電車で戻って歩き出せたことであった. 追分～下諏訪, 大井～太田のような山道区間を除けば行程はかなりフレキシブルであった.

自治体により中山道の周知は様々であった. 「歴史の道 中山道」の道標が立っているところもあれば, 全く無関心な町もあった. 中山道の保存にもっと尽力して欲しいものである. 中山道を歩く人がさらに増えるように.

3シリンダ機関車

前回書いた3シリンダ機関車の中央シリンダの制御シャフトの続きで, リンク機構はベクトルで考えると明瞭なことが分かった.

120度づつ離れたベクトル a, b, cがある. aに2 to 1 leverを接続すると, 相手の端は逆向きで大きさが半分のベクトル a'になる. これはちょうどbとcの中間なので, a'を中心にしたbの動きはちょうとcになるわけだ.

3シリンダ機関車

時間があったので鉄道博物館へいった. いつも何か発見があり楽しい.

国産の3シリンダ機関車C53の走り装置が置いてあった. 3シリンダ機関車は, 両側のシリンダは普通の蒸気機関車のとそっくりだが, 外から見えない台車の中央に3つめのシリンダーがある.

その中央のシリンダからはクランク軸になっている動輪に力を伝達する. 問題は中央の制御シャフトがどうなっているかであるが, 実に意外な構造であった. 両側の滑り弁を動かす制御シャフトが機関車の先頭の方へ突き抜け, それをリンクでつないで中央の制御シャフトを動かすのであった.

googleで探したら,
http://www.watercressline.co.uk/tw/pics/bitn2to1.jpg
にその図解があった.

こういう計算をしてみた. 一方の制御シャフトの動きをa, 他方のそれをb, 中央のそれをcとする.

a=sin x
b=sin (x+2π/3) = sin x cos 2π/3 + cos x sin 2π/3=-1/2 sinx + √3/2 cos x
c=sin (x+4π/3) = sin x cos 4π/3 + cos x sin 4π/3=-1/2 sinx - √3/2 cos x

これを眺めると a + b = -c とわかる.

たしかにベンツのマークのような, 120度ずつはなれたベクトルをa, b, cとすると, a+bはcと反対向きのベクトルになるから, 納得できる.

さて上の図解によると, 下のvalve spindle linkの動きをfulcrum pinを中心にして2 to 1 leverで上へ伝達する. その点を中心として, 上のvalve spindle linkの動きをequal leverで中央のvalve spindle linkへ伝えている.

従って, 下の動きをaとすると, 2 to 1 leverの上の動きは -1/2 aになる. 上の動きをb, 中央の動きを x とすると, bとxの平均が-1/2 aなのだから (x + b)/2 = -1/2 a. したがって x = - (a + b)となり, ちょうどcの動きとなっている.

すばらしい.

銀河

京都へ出張した帰路, 3月14日で運転をやめる寝台急行銀河を利用した. 以前銀河で下ったときの乗客は半分以下で, 車掌といつまで運転するのかなぁと話したりした. 運転打ち切りの直前ともあって今回は満員であった.

寝ながら聞く鉄橋を渡る響き, すれ違う列車の音, まれに停車する駅の物憂いアナウンス, 「おはようございます. 何月何日何曜何時何分です」という車掌の声で朝を知る. 独特の寝台列車の朝の雰囲気である.

子供のころから寝付かれない夜は今寝台列車で旅していると考えるとすぐ眠れた. 「銀河鉄道の夜」も何度も読んだ.

最近利用する寝台列車は函館からの帰りの北斗星ばかりだが, 利用回数の多いのは銀河であった.

新幹線を併走する路線で寝台列車を利用するのは, 物好きといわれても仕方ない. 能登は残っているし, これからは金沢へ向かう能登と函館から戻る北斗星をときどき利用するだけになりそうだ.

ついでに
英字新聞を読んだ人からアメリカでは洪水で何万という人が寝ていて流されたと聞いたがこのsleeperは寝ている人ではないらしい. sleeperは寝台車であると同時に線路の枕木のことでもある.

knuthのGPS

電通大での講義の時に，GPS(general problem solver)の話をした．Algol 60の名前呼びの機能を最大限活用するものだ．


real procedure gps(i,n,z,v); real i,n,z,v;
  begin for i:=1 step 1 until n do z:=v; gps:=1 end


つまり4つの引数i,n,z,vはvalueと宣言していないからすべて名前呼びである．

これを使うプログラムは

i:=gps(i, if i=0 then -1.0 else i,p,
　if i=1 the 1.0 else
　　if gps(a,i,z,if a=1 then 1.0 else
　　　entier(a)×(entier(i)÷(entier(a))=entier(i)∧a<i then 0.0
　　　　else z)=z then
　　　　　(if p<m then=p+1 else i×gps(a,1.0,i,0.0)) else p)

プログラムを実行するとpにm番目の素数が入る

最初の行はiを制御変数にしてfor文をまわす．iは1になっているので終点はiということ，すなわちiはどこまでも増え続ける．forループで実行するのは変数pへの代入で，まずはiが1なので1.0が入る．
次はiが2になって代入文へくる．今回からpへ代入する値はelse以降3行目からになる．
そこでgpsを呼ぶ．今回は制御変数aで1からiまで増える．代入先はzで初回は1.0が入る．その後aは1ずつiまで増える．代入される値は4行目で計算する．ここでentierはAlgol 60の基本関数で小数点以下を切り捨てる. ÷の記号は商の小数点以下を切り捨てる除算.　したがってここではiがaで割り切れるかを見ている.　割り切れてしかもaがiより小さければiには約数がある，素数ではない.　そうならzにoを入れ，そうでなければzのまま. zには最初1を入れたが途中1つでも約数があると最後は0になっている.
一方このgpsから脱出してきたときの値はgpsの宣言の最後にあるように1である. zがそれに等しいというのはiが素数であったことである. そこで最後の行のelseより前を実行してその値をpに代入する. pがmより小さければp+1だからpはmになるまで素数のたびに1ずつ増える. 素数でなければp, つまり不変.
ということでm番目の素数が見つかったところで中ほどのelseの次に来る. iには今m番目の素数が入っている. gpsの値1を掛けてpに入れる. gpsの仕事はiに0を代入して外側のgpsのループを停止させることである. 問題は最後の乗算でiを先に評価すればよいが, gpsを先に評価すると大切なiの値が変わっていることだ.
まだ他にも微妙な点はあるが，いちおうこのようなプログラムだ.

計算機プログラムの構造と解釈

この冬学期，電気通信大学の大学院で記号処理特論という講義をしました．この学期には「計算機プログラムの構造と解釈」の4章，5章を説明しました．90分の講義15回でカバーするのは結構大変です．

私は実装に興味があるので，使い方の例題はほとんど無視してしまいました．また問題にも触れませんでした．

相当なスピードで進んだので，消化不良の学生さんもいたと思われます．

希望的にはＭＩＴのサイトからプログラムをダウンロードし，さまざまな例題でプログラムを走らせ，体験を通して理解して欲しかったと思います．（そうした人がいなかったとは言い切れません）

私個人は4章，5章が1350分で講義できたという経験が貴重です．講義の機会をいただいた電気通信大学に感謝しています．