2015年12月11日金曜日

菱形六十面体を織る

前回のこの題名のブログの終の方に「ちょっとややこしいのが糊しろの部分で, 正六面体のように簡単ではない. その説明はまたの機会にしたい.」と書いたその部分の説明だ.

下の絵を見てほしい. 線がごちゃごちゃ描いてあるが, 前回の最後の図の各帯が菱形六十面体のどの面を形成しているかを示している.



例えばAの帯は, 真中少し下のAの回りのオレンジ色の線に対応している. Aの, 時計でいえば1時の方向に0の面があり, そこから右下へ進んでαの辺から外へ出, 少し下のαの辺から展開図の戻り, 一つおきに1,2,3の面を通り, θから再び外へ出, 上のθの4から戻ってくることを示している.

緑で示すBの帯は, 5,6,7,8,9だが, 5はIの五角形の下の方にあり, 右上に進み, γから出て, 右下のγから戻り, Bを取り巻くように6,7,8,9が見える. Bの場合, 外へ出るのは1回である.

KとLは, 外へ出ないから遥かに簡単である.

ところでこれらの閉曲線には, その色で塗りつぶした三角形と, 白抜きの三角形が書き込んである. 実はそこが糊しろである. 白抜きの方が出発点. 線上にある三角の頂点の方へ進む. これが前回の帯の絵の上端に相当する. Aで言えば0の上の白い部分(下の49と重なる部分)であり, 塗りつぶしの方が下端になる. 線上の三角の頂点の方から下端に到る.

これを眺めると, K, Lは別として, 他の帯は上の五角と下の五角を貫いていて, A,B,C,D,Eは下では円弧, 上では直線であり, F,G,H,I,Jは逆になっている. K以外の糊しろはすべて上の五角の方にある.

つまり下の五角は出来上がった正十二面体の下半分の鉢を構成し, 上の五角は上半分を構成していたわけだ.

底の方から組み立てて, 糊しろが込み合う上面で組み立てが完了するような構造になっていた.

だから組み立て方は, A,B,C,D,Eをほぼ中心で重ね合わせ, その周囲をKの帯で固定し, その辺からF,G,H,I,Jも加わって下の鉢を作り, 中段を過ぎて上の周囲の面に至り, それをLで補強した後, F,G,H,I,Jを上の面に相互に差し込むのである. (糊しろといっても糊をつけて貼るのではない.)

元々の帯が巧みに設計されていて, この図を描くにも, 円弧と円と直線の基本の曲線群を用意してしまえば, たちまちにして完成した次第だ.

2015年12月3日木曜日

EDSACのプログラム技法

EDSACのプログラムを集積したウェブページがある. (詳しくはそこの "To browse the collection of software, click this link."をクリックする.) ここには私の書いたエラトステネスの篩も置いてあって嬉しい(21行目). そこにDavid J. Wheelerの素数のプログラムがあった(15行目). Martin Campbell-Kellyの書いたEDSACの文献(Programming the EDSAC: Early Programming Activity at the University of Cambridge, IEEE Annals of the History of Computing, Vol 1, No. 2.)によると, EDSACには「EDSACの本」で有名なイニシアルオーダーの前にもう一つのイニシアルオーダーがあって, 1949年5月6日の運転開始の日にはその方が使われていた. 同年6月22日にケンブリッジで開催された計算機会議でEDSACが公開され, その時, Wilkesの書いた平方の表とWheelerの書いた素数のプログラムのデモが行われた.

Wheelerは1927年生まれだから, 22歳のWheelerが当時どういうプログラムを書いていたかには非常に興味がある. これはもう読むしかない.

ところがその内容たるや
[Primes]
T107SO92SO93SO94SS5ST6SO95SO95ST7S
A96SR4SS97SE42SL4SA97SG45SS98SG100ST7S
A97SA4ST97SH97SN97SL64SL64SA96SE39ST7S
A96SU1SA4ST96SA99ST97SS88ST7SH91SA1S
E72SV91SS89SE71SA89STLOSH90SV1SL4S
TLA7SA98SG67SA6SA4SG36SE33SP2SP500S
JSP16S#S@S&S!SP2LP1LPLP1L
T7SA99ST97SA4SA96ST96SE39S
だけなのだ. これを読む手がかりはここにある.

読んでみると, やはり面白いプログラムであった. 特にそのプリントルーチンには感心した. 以下にその解説をしたい.





31番地 T107S. このプログラムは最初のイニシアルオーダーInitial Order1で読み込むが, 読み込まれるプログラムが106番地まであるなら, 最初をT107Sとするのであった. つまりイニシアルオーダーには格納命令があり, 最初はT31Sになっていて, 31番地から命令を格納する. つまり, T31S, T32S, ..と変わっていく. そして格納作業の終りは, 格納命令から31番地にある命令を引き, それが正であるかを見る. T107Sを引くと0(正)になるので31番地へ進んでくる. このT命令はアキュムレータのクリア用にも見えるが, すでに0になっている.

32番地は92番地にある#S(figure shift)を出力する. 以下のプログラムが数字を出力するので, その準備.

33,34番地では93, 94番地にある@S(cr), $S(lf)を出力. 次のS5Sはイニシアルオーダーの5番地にあるP5S(10*2^-16)を引き, -10を作り, 6番地(m[6])に入れる. 6番地は -10,-8,-6,..と変わり, 1行に5個の素数を印刷する.

37, 38番地. !S(space)を2個出力. (このプログラムのあちこちにある)T7Sは(67番地のものを除いて)アキュムレータをクリアする.

このプログラムでは素数の候補pから3,5,...と奇数(d)を引き(以後 dを減数という), 0になるかどうかを見ている. (前にも書いたようにEDSACには除算命令はない.) pは96番地, dは97番地にあって, それぞれ5と3に初期化されている.

40から46番地はそのpからdを引く. ただ初めは16*dを負になるまで引き, 次にdを正になるまで足す. これでpをdで割った剰余が得られる. 実際には16*dを引くのではなく, pを4ビット右シフトして引き, 4ビット左シフトして足す. (おおざっぱにいうと単にdを引くより, 16倍加速されている. こういう剰余の取り方はEratosthenesの篩でも使える.)

47番地. EDSACには零ジャンプがないので, 正になったところで1を引き, 負なら零であったということで, 合成数の判定をして100番地へ進む.

剰余が0でなければ, 50,51,52番地でdを2増やし, 53番地からH97Sで乗算レジスタにdを置き, N97S, L64S, L64Sで-d^2を作る. それにpを足し, 減数の二乗が素数候補より小さければ39番地へ戻り, 次のテストに掛かる.

減数が大きくなって素数と判定出来ると, 60番地から先で素数の出力にとり掛かる.

U1Sで素数を1番地に入れ, A4Sでpを2増やして63番地で96番地に入れておく.

素数は4桁で出力するつもりで, 66番地で-4を7番地のカウンタへ入れる.

68番地から76番地は1桁の出力で, 1番地に1234があるとして説明すると, 68番地H91Sで乗算レジスタに2^-11を入れる(32*2^-16=2^-11).

アキュムレータに1234を足す. 72番地へ飛んで89番地にある1000を引く. アキュムレータは234(*2^-16)になり, 正なので71番地へ行く. V91Sだからアキュムレータは234*2^-16+2^-22になる.

再び1000を引くと今度は負になり, 74番地で1000を足し, アキュムレータは234*2^-16+2^-22に戻る.

これを0番地の長語に入れる. つまり左半分の1番地は234, 右半分の0番地はファンクション部に1となる. そこでo命令は0番地を出力するから1が印字される(EDSACの文字コードは2012年9月28日のこのブログ参照).

77番地から80番地は234を10倍する. H90Sで乗算レジスタに10*2^-4を置き, V1Sでアキュムレータは234*10*2^-4になり, 次の左シフトで2340になり, 1番地に戻される.

81番地からは7番地のカウンタを1足し, まだ4桁出力していなければ, 67番地へ戻る.

次は2340だから, 71番地のループは2回実行され, 右半分の方は2*2^-11になって2が印字され, 1番地は3400になる.

4桁出力すると84番地へ来て, 6番地のカウンタにA4Sで2を足し, 5個未満なら36番地, 5個終わっていたら33番地へ戻る.

一方合成数の時は100番地へ進み, 101, 102番地でdを3に戻し, 103,104,105番地でpを2増やして39へ戻る.

この解釈が正しければ, 期待される出力は
  0005  0007  0011  0013  0019
  0019  0023  0029  0031  0037
...
なる筈であった. サイボウズ・ラボの西尾君がEDSACのシミュレータを持っているというので, 試してもらったらその通りであったとメイルが来た.

このプログラムの疑問は素数候補が5から始まることである.

減数を3から始めると素数候補の3は割り切れるから, 3は素数にならない. 3を素数にするには2から割り始めなければならないからなのだ.

2015年11月30日月曜日

菱形六十面体を織る

ネット上に山口大学理学部数楽工作倶楽部の菱形60面体編みというページがあった.

なんだか面白そうだとさっそく作って(編んで)みた. 下の写真がその完成品.



編むというのは糊を使わず, 紙の摩擦だけで固定されているからである. そんなことで菱形六十面体が作れるのかとうのが興味の中心である.

面を1/10にして正六面体で考えてみると, 下の図のように正方形を5枚繋げた(ペントミノの`I'の形)を環状にした帯を3個作り, 左から3個のような方向に置き, 左からx, y, zということにする. xとyの面が重なる手前左方向ではyがxの外に, yとzが重なる右方向ではzがyの外に, zとxが重なる上方向ではxがzの外になるように組み立てると, 右端のような六面体が出来る. いわゆるみすくみ状態になって下の帯を上から抑える構造に出来ている.



菱形六十面体もこの方針で作られいるようだ.

私は元の型紙をダウンロードしたのではなく, PostScriptで書き直し, 少し厚手の紙に印刷し, 山折り谷折りをした後, カッターで帯を1枚ずつ切離し, 裏の番号を写してから組み立てに取り掛かった. 小さいクリップを5個使い, 組み合わさった場所を仮止めしながら進む.

意外な方向に帯が折れていくのに戸惑いながら完成したのがこのブログの最初の写真である. その後, 構造を理解しようと思い, 各面に番号を振ったり分解したり組みなおしたり, 結局3回作った. もちろん後ほど早く完成する.

さて菱形六十面体は, 私の5年ほど前のブログ にあるように, 正十二面体が基本になっている. その12個の正五角形の各面を5枚の黄金比の菱形に置き換えたもので, 面も稜も凹んでいる. 正十二面体の各面を(変形した)菱形にした展開図が下だ. 12の面には白抜きでAからLまで記号があるが, これは元の型紙の12本の帯の左から右へA,B,C,...としたものに対応する.



真ん中辺のAの正五角形の外周に, オレンジ色の0,1,2,3,4が見えるが, これが型紙の左端の帯で, 0と1の間は9の下を潜り, 1と2の間は50の下を潜り, 2と3の間は22の下を潜り,...という風になって, 0はJの面の, 1はBの面の, 2はKの面の一部を構成する. つまりAの帯はAに接続する5個の面を作るように出来ている.

そうしてみるとAの帯は, 0, 9, 1, 50, 2, 22, 3, 28, 4, 49の菱形を外に出たり下に潜ったりしてして立体を作ることが分かる. というより, 型紙の各帯の外に面する菱形に, Aには上から0から4, Bには5から9, Cには10から14,.. Lには55から59と番号を付け, その番号を展開図に書き込んだものである. (外周のギリシア小文字は対応する辺同士を示す.)

この番号を型紙に写したのが下の図だ.



A,B,C,D,Eの帯の上から4段目が左から50,51,52,53,54. F,G,H,I,Jの帯の最下段が左から55,56,57,58,59になっているのは, それぞれKとLの帯の下を潜って, KとLの面の周囲を固めているのである.

ちょっとややこしいのが糊しろの部分で, 正六面体のように簡単ではない. その説明はまたの機会にしたい.

もしかするとこの編み方はこちらが先かもしれない.

2015年11月28日土曜日

Zuseの計算機

Zuseの計算機Z1は金属板の移動で論理素子を実現していて, その金属板は重なっていてるから, これで上手く動くか心配になる.

現に, 戦後再生されたZ1はZuseの存命中は稼働していたが, その死後にはもう動かなくなっているらしい.

さて, 前回の基本素子に続けてAND, OR, XORの回路は次のようだ. 前回のリレーもゲートだからANDでもあるが, 今回のは2つの入力のANDでゲートを制御するものである. (OR, XORも同様)

先ずはAND. 下の2個は入力の制御板, 上の右が出力の受動板, 上の左がクロック信号に相当する能動板である. 然しANDにはもう1個, 中段に遊動板というものが存在する. この図は入力A, Bが共に0の状態で, 棒は2本とも下の位置にある.



この時能動板が右に動いても, 能動板の棒の左に隙間があるので, 棒は動かず, 受動板も動かない. つまり出力も0である.

A=1, B=0の時は, Aの棒が上にいるので, 能動板に従って棒は右に動く. すると遊動板も棒に押されて右に動く. しかし右側の棒はB=0なので下にいて, 遊動板が右に動いても棒が右に動くことはなく, 受動板も動かない.

反対にA=0, B=1の時は遊動板の右の棒は上の位置にあるが, そもそも遊動板が動かないので, 受動板も動かない.

両方が1なら遊動板も動き, それに押されて右の棒も右に動き, 受動板も動いて出力が1になる.

ORは遊動板もないのでもっと簡単だ.



下の2個が制御板である. 上の左が能動板. 右が受動板である.

入力が共に0, 2本の棒が下にあると, 能動板が右に動いてもどちらの棒もそのままなので, 受動板は動かず, 出力も0である.

いずれかの入力が1なら(共に1でも), 少くても1本の棒が上にあがり, 能動板が右に動くと上にあがっている棒に押されて受動板が右に動く.

排他的論理和XORは一番面倒である.



下の制御板, 上左の能動板, 上右の受動板の他に中段に遊動板が2枚ある.

まずA=0, B=0の時. この図の構造で能動板が右に動く. すると上の遊動板は右に動くが, 右の棒の左に隙間があるので, 右の棒は動かない. 下の遊動板は左の棒の右に隙間があるので, これは動かず, 右の棒も動かない.

A=1, B=0の時に能動板が動くと, 上の遊動板は動かないが, 下の遊動板は動き, 下の遊動板の右の棒も動いて受動板が動く.

A=0, B=1の時の能動板が動くと, 上の遊動板が右に動き, それに従ってその右の棒も動くから受動板も動く.

A=1, B=1の時は, 上の遊動板は動かず, 下の遊動板は動くが, 右の棒の左の隙間のせいで右の棒は動かず, 受動板は動かない.

なかなかよく出来ている.

2015年10月27日火曜日

Zuseの計算機

Konrad Zuseはベルリン工科大学で, 電気工学や機械工学でなく土木を学んだが, 1930年代の学生時代にすでに計算機を作ることを考えていたといわれる.

Z1からZ4までのZuseの計算機は全く独特なもので, 特に1936年〜1938年に作られた全機械式のZ1に私は興味がある. すべて金属の板, 金属の棒のたぐいでできた浮動小数点の加算器, 制御装置, 記憶装置はまったくユニークであり, その構造を調べたくなるのは人情であろう. 1936年といえば, あのTuringの論文の出た年である.

しかしこのZ1は大戦中に連合軍の空襲で破壊され, Zuseが戦後に再構築したZ1が存在し, その写真をみることができる.

再構築の話はここに詳しい.

今回からしばらくZuseの計算機についてブログを書きたい.

アメリカのComputer History MuseumにZ1の一部の写真がある. その基本素子の構造は次のようだ.



3枚の金属製の板に穴が開いていて, その穴に棒が通っている. 金属板はこの図では下から制御板, 受動板, 能動板という名称だが, 板の重ねる順は別に問題ではない. 問題は板の動く方向である.

制御板はyと書いた奥方向に動き, 手前方向に戻る. 受動板と能動板はxと書いた右方向に動き, 左方向に戻る. それぞれ現在の位置が0の位置で, 動いた先が1の位置である.

制御板が図の位置にあるとき, 能動板が動くと, 棒の左側に穴が続いているので, 棒は押されることはなく, 従って受動板も動かない(0のまま).

制御板が奥方向に動いたとすると, 棒が奥方向へ押され, 受動板と能動板の穴の奥の方に移動する. ここで能動板が右に動くと, その辺の穴は左に余裕がないから, 棒も右に押され, 従って受動板も棒におされて右に動き, 1の位置に来る.

つまり制御板が0なら受動板も0, 制御板が1なら受動板も1という情報伝達ができる.

右に描いたのは電気的なリレーで, 下のコイルに電流が流れると, 上の回路が閉じ, 回路に電流が流れる. コイルの電流が切れると回路の電流も切れるのと同じなので, この基本素子はリレーだという人もいる. 電気的リレーは, 応答が瞬時だが, Zuseの素子では, 制御板より能動板が遅れて動くので, 出力の受動板の反応も遅れる.

さらに入力方向と出力方向は90度ずれていることにも注意しなければならない.

ところで受動板が出力を出したら能動板は戻ってよいかといえば, それは駄目で, 能動板が戻ると受動板も戻ってしまい, これは次の論理の入力になっているからまだしばらく止めておかなければならない.

入力の制御板はすでに自由である. という次第で, 入力が0(左側), 1(右側)で信号を伝える分解図が次である.

この図では下から制御板, 能動板, 受動板の順になっている.



Cycle 0. 初期状態 受動板も能動板も定位置(0の位置)にいる. 着目するのは制御板と能動板の間が広く, 制御板と能動板の縦の縁がずれていることだ.

Cycle 1. 右側は制御板が上って1を入力した. 左はそのまま. 右の図の棒は上側へ移動する. 制御板と能動板の間が狭く, 両板の縦の線はずれている.

Cycle 2. まん中の能動板が右へ動く. 右の図の棒は左へ移動し, 受動板も右へ移動する(1を出力) 制御板と能動板の間が狭く, 両板の縦の線は揃っている.

Cycle 3. 制御板が下へ戻る. 制御板と能動板の間が広く, 両板の縦の線は揃っている.

Cycle 0. 能動板が左へ戻り, 初期状態になる.

つまり, 制御板が動き, 能動板が動き, 制御板が戻り, 能動板が戻る という動作をくりかえす.

これで90度曲ったが, 1で入力した信号は2で出る, と同時に次の入力になるから, この4 clock の間に4回の論理演算が出来たことになる.



上のcycle 1の右と同様な絵が左上に見える. Cycle 1の特徴が見える. 下から来た縦長の制御板が上に移動して, 1を入力したところだ.

右上の素子は横向だがCycle 0の所である. 右下は制御板の能動板の間に隙間があるから, 制御板が戻った所でCycle 3である.

左下はCycle 2で能動板が動き, 受動板を動かし, 出力と同時に右上の素子に入力した.

この動画がここにある.

Zuseのこの素子は4クロックでフリップフロップを実現しているが, 60年程前にあったパラメトロンは3拍励振といって, I相の論理素子の出力をII相へ入れ, II相の出力をIII相に入れ, III相の出力をI相へ戻すという3クロックの論理回路を構成していた. 久し振りにそんなことを思い出した.

2015年9月25日金曜日

ARMのプログラム技法

ARMには除算命令がないので, 前のブログCalの時はdivmodという簡単な除算サブルーチンを用意した. 通常はライブラリの除算ルーチンが組込まれるらしいが, 私としては除算サブルーチンを書くのが趣味の内だ.

2008年10月28日のブログに書いた「個人用電卓」の除算もそうだが, 私は剰余を除数と符号を合せる除算が好きなので, 今回も当然そういう設計にした.

400を15で割ると, 商は26, 剰余は10になる. では-15で割るときはどうするか. 商は-27, 剰余は-5とするのである. 通常割切れるときは剰余は0だが, これは正なので, 負数, 例えば-15で割ったときにはそういう剰余はでない. 剰余は-15になる.

400割る-20は商が-21, 剰余が-20とする. -21× -20=420, それに剰余が-20だから, 被除数は400というわけだ.

引き放し除算の方法はこうだ. ARMのレジスタは32ビットだが, 4ビットだと思って書いた図が下だ.

左は55割る7, 右は47割る7.



行0. 被除数が二進法で置いてある. それを左へ1ビットシフトする.
行1. 被除数と行2にある除数の符号を較べる.
行2. 符号が同じなら被除数から除数を引き, 異なれば足す.
行3〜12. 被除数を左へ1ビットシフトする.
シフトする前の被除数と除数の符号が同じなら被除数から除数を引き, 被除数の最下位に1を置く. 異なれば足す.
行13. 右のレジスタを2倍して1を足す. 左のレジスタに剰余. 右のレジスタに商がえられる. 剰余と除数の符号が異なれば, 剰余に除数を足し, 商から1を引く.

これと同じことをARMで実行したのが下のプログラムである.

プログラム0


行00 .data ここからデータ領域という指示
行02〜04 テストデータ
行05 printfの書式
行08 ここから除算サブルーチン. 被除数は上位がr1, 下位がr0, 除数がr2にある.
行08,09 r1,r0を繋げて左シフト. r0+r0をr0に置き, addsでキャリーをcpsrに 置く. adcでそのキャリーを足しながらr1+r1をr3に置く. adcsなので オーバーフローがあればセットする.
行10 被除数の左端の2ビットが01, 10だとオーバーフローが起き, bvsでジャンプ.
行11 除数と被除数の符号を較べる.
行12,13 符号が同じなら被除数から除数を引く, 異なれば足す.
行14 その加減算の前後の被除数の符号を較べる. 変っていなければオーバーフロー.
行16 カウンタr4に30をセット.
行17,18 被除数を1ビット左シフト.
行19 符号を較べる.
行20〜22 同じなら除数を引き, r0に1を足す. 異なれば足す.
行23,24 カウンタを減らす.
行25,26 補正
行27〜29 補正
行30 サブルーチンから帰る.
行31〜33 オーバーフローのとき, 除数が正なら-1, 負なら0を持って帰る.

プログラム1


これはドライバのプログラム.

クルティカルな例を実行したのがこれだ.

被除数            除数     商       剰余
3fffffff 7fffffff 7fffffff 7fffffff 7ffffffe
c0000000 80000000 7fffffff 80000000        0
c0000000 7fffffff 80000000 7fffffff ffffffff
3fffffff 80000000 80000000 80000000 80000000
Schemeで検算してみる.
(number->string
(+ (* #x7fffffff #x7fffffff) #x7ffffffe) 16)
=> 3fffffff7fffffff

(number->string
(* #x7fffffff #x-80000000) 16)
=> -3fffffff80000000 (= c000000080000000)

(number->string
(+ (* #x-80000000 #x7fffffff) -1) 16)
=> -3fffffff80000001 (= c00000007fffffff)

(number->string
(+ (* #x-80000000 #x-80000000) #x-80000000) 16)
=> 3fffffff80000000
うまくいっているようだ.

2015年9月16日水曜日

ARMのプログラム技法

数日前のブログではARMのアセンブリ言語でカレンダーのプログラムcalを書いてみた. 今回はおなじラズパイで動くBCPLで同様のプログラムを書こう. (ARMのプログラム技法という題は多少おかしいが.)

BCPLはMartin Richardsが1967年頃に設計した言語である. その元はケンブリッジ大学とロンドン大学がコンパイラ記述用に開発したCPLで, Combined Programming Languageの略である. そのBasic版がBCPLということになっている. またBCPLからC言語が生まれたこともよく知られている.

BCPLをケンブリッジで開発した後, Martin RichardsはMITのProject MACに滞在したので, MITでもBCPLは使えるようになった. 私がBCPLに出逢ったのは1973年から1年MITにいた時であった.

BCPLは最初の設計から50年近く経ったが, まだ保守されているのに驚く.

下のcal.bはBCPLでcalを書いたリストである.



ARMのcalとほとんど同じに出来ているので, それを読む注釈と思って眺めて欲しい. 00行目 C言語のinclude のようなものだ.
02行目 C言語のmain()のようなものだ. ただVALOF { ... RESULTIS ..} の形に なっている.
03行目 変数yearとmonthを定義. BCPLには型がない. 初期値もこう書く流儀らしい.
04行目 他の変数も定義
05,06行目 配列の定義 VEC 11で0から11まで12個の場所が取れる. 初期値の設定は 別にやるらしい. 14行目から18行目で設定した.
08行目から13行目 yearとmonthを引数から取り込むにはこう書くらしい.
19行目 Gregorio暦では週日は400年で繰り返すから, 西暦を400で割った剰余r400だけ を考えることにする.
21行目まで r400を100, 4で割った商と剰余を求める.
22行目 r100≠0をBCPLではr100¬=0と書いたりするが, ASCIIには¬がなく, ~=0でもいいようだ.
23行目 閏年ならrsの1月と2月, msの2月の値を修正する.
24行目 aはその月の1日の週日. nは1週間を数える変数である.
25行目 1日のある週の最初から6週目の最後までFOR分で繰り返す.
26〜28行目 1日から月末までの日なら出力する. それ以外なら空白を出力.
29,30行目 nを下げ0になったら改行する.
31行目 0を持って帰る.

これを実行するには, raspberrypi:~$でcintsysに入る

0.030> bcpl cal.b to cal
bcpl cal.b to cal

BCPL (10 Oct 2014) with simple floating point
Code size =   436 bytes of 32-bit little ender Cintcode
0.110> cal 2015 9
cal 2015 9
        1  2  3  4  5
  6  7  8  9 10 11 12
 13 14 15 16 17 18 19
 20 21 22 23 24 25 26
 27 28 29 30

0.020> cal 2015 10
cal 2015 10
              1  2  3
  4  5  6  7  8  9 10
 11 12 13 14 15 16 17
 18 19 20 21 22 23 24
 25 26 27 28 29 30 31

0.020> 
とまぁこんな具合に動く.

2015年9月10日木曜日

ARMのプログラム技法

前回のブログでCalのプログラムの核心の部分がうまく書けることが分かったので, 全体を書いてみた. やや長いが興味があれば追ってみて欲しい.

プログラム0


行00 .data ここからデータ領域という指示
行01,02 年と月を格納する場所 4バイトずつ
行03 1月から順に次の年初の曜日とその月の1日の曜日の差 12月は31日なので 4週と3日余るが, -3mod7=4のような値のリスト. 以下rsという
行04 各月の長さ. 以下msという
行05〜09 scanf, printfの書式
行10 プログラム部分の開始
行11〜17 除算サブルーチン r0をr1で割り商をr0, 剰余をr1に置く

プログラム1


行19 主ルーチンの入り口, 帰り番地をr8へ退避する
行20 yearを入れる場所をr1へ
行21 monthを入れる場所をr2へ
行22 scanfの書式の番地をr0へ
行23 scanfを呼ぶ
行24 データ領域のベースアドレスをr3へ
行25 yearをr0へ
行26 400をr1へ
行27 除算
行28 400で割った剰余R400をr4へ
以下 R400を100と4で割った商と剰余を求める
行38 100で割った剰余を0と比較
それが0か否かにより, 400と4との剰余をr1に置く(条件つきmov)
行41 r1が0なら閏年(閏年の判定には100で割った剰余R100が0なら400で割った剰余R400が0, 0でないなら4で割った剰余R4が0かで判定できる)

プログラム2


行42〜47 r1=0なら閏年なので, rsの1月2月とmsの2月を修正する
行48〜58 その月の1日の週日(a)は(1+R400-Q100+Q4+rs[月])mod 7
行59 前回のブログのr5の値, 1-a
行60 前回のブログの#6に相当する値, 43-a
行62 前回のブログの#3に相当する値, その月の日数

プログラム3


出力用の値が設定できたので出力に移る.
行63,64 曜日の名を出力
行65 r4は1週間を数えるカウンタ
行66 日付が1より大きいか
行67 そうなら最後の日付より小さいか
行68 日付をr1に置く
行69,70 比較の結果でいづれかの書式をr0に置く
行72 日付を1増やす
行73 週のカウンタを減らす
行74〜76 1週間が終わったら改行しカウンタをリセット
行77,78 最後になっていなければl0へ戻る
行79 lrを戻す
行80,81 プログラムの戻り値を入れて終わる
行82〜87 定数の番地

このプログラムでは, データの場所にはラベルがあるが, プログラムの部分には殆どラベルが見られない. divmodはサブルーチンの入り口, mainはプログラム全体の入り口, あとカレンダー出力のループのl0があるだけである. 従ってジャンプ命令もほとんどない. そういうことではARMのアーキテクチャは気持ちがいいといえる.

このプログラムをRaspberry Piで動かすには,
as -o cal.o cal.s
gcc -o cal cal.o
./cal
で起動し
2015 9
のように年と月を入力すると
./cal
2015 9
 Su Mo Tu We Th Fr Sa
        1  2  3  4  5
  6  7  8  9 10 11 12
 13 14 15 16 17 18 19
 20 21 22 23 24 25 26
 27 28 29 30

とカレンダーが得られる.

このプログラムはhttp://www.iijlab.net/~ew/cal.sに置いてある.

2015年9月8日火曜日

ARMのプログラム技法

1ヶ月ほど前のブログで「数式お絵書き」を書いたのは, 手元のRaspberry PiでMathematicaが動いていたからだ.

Mathematicaも面白いが, Raspberry PiはそのCPUがARMなので, アセンブリ言語でプログラムを書くと, ARMのアーキテクチャがよく分るではないかと, 最近ではもっぱらアセンブリ言語のプログラムを書いている.

「ARM」という単語を最初に聞いたのは, 1995年頃の英国ケンブリッジでであったと思う. ケンブリッジ大学の計算機研究所のWilkes先生を訪ねたのだが, もう大学にもオリベッティの研究所にも出勤されてはいなかったので, オリベッティ研究所の所長のAndy Hopperさんに会って雑談している時, ARMといういいCPUがあるよと言われた. (この訪問のことは, bit 1996年1月号のaleph zeroに書いた.)

昨年8月, 偶然みつけたウェブページのことをツィッターに書いた.
50年ほど前にBCPLを開発したMartin Richardsが「青少年のためのRaspberry Pi上のBCPLプログラミング入門」を書いている(http://goo.gl/rXOKy2 ). 1つの言語に半世紀も情熱を持ち続けるとは見上げたもの. そのうち読んでみたい.
(「青少年のための...入門」と書いたのは, 元の題が「Young Persons Guide to BCPL Programming on the Raspberry Pi」であり, Benjamin Brittenの曲「The Young Person's Guide to the Orchestra」が日本では「青少年のための管弦楽入門」といわれているからだ.)

このツィートの後, 研究所でRaspberry Piを何個か購入し, 使えるようになってきた. またARMそのものにも興味も出てきた.



ARMでプログラムを書くに際して, 最低限の知識は次のようだ.

ArmはRiscだが, 命令幅が32ビットで, 水平型マイクロプログラムの様相を持つ. 16ビット命令のモードもある.

記憶装置とのデータ転送はldr, strだけ.

CMPがCPSR(Current Program Status Register, NZCVビット)をセットするほか, 加減乗算命令でも, addsのようにsを付けてCPSRをセットできる.

命令には, eq, ne, pl, miなど条件を付け, 実行するしないが指定できる.

このようなことを念頭に置いてプログラムを書く. その前後に指定された個数の空白を置き, aからbまでの数字を順に出力するにはどうするか. これはカレンダーのプログラムを書くときに必要になる. 例えば今年の9月のカレンダーは, cal 9 2015 で見ると
% cal 9 2015
   September 2015
Su Mo Tu We Th Fr Sa
       1  2  3  4  5
 6  7  8  9 10 11 12
13 14 15 16 17 18 19
20 21 22 23 24 25 26
27 28 29 30

なら, 先頭に2個の空白と最後に3個の空白を置き, 1から30まで印字するプログラムを書けばよい. いきなりプログラムではなく, テストプログラムとしては, 例えば前に3個, 後に2個の空白を置き, 1から3まで出力するプログラムは普通ならこうなろう.
(do ((i 0 (+ i 1))) ((= i 3)) (display "   "))
(do ((i 1 (+ i 1))) ((> i 3)) (display i))
(do ((i 0 (+ i 1))) ((= i 2)) (display "   "))
ARMの機能を活用するとこうなる.
.data
        .balign 4      /*4バイト境界を合わせる*/
a:      .word -2       /*出力する整数の初期値 2*/
ret:    .word 0        /*帰り番地の退避場所*/
s0:     .asciz ". *"   /*printfの書式*/
s1:     .asciz ".%2d"
s2:     .asciz "\n"
.text
        .balign 4
.global main
main:   ldr r0,reta    /*帰り番地を一時退避*/
        str lr,[r0]
        ldr r0,aa
        ldr r5,[r0]    /*出力する整数をr5へ*/
l0:     cmp r5,#1      /*出力する最初の数a*/
        rsbgts r0,r5,#3/*出力する最後の数b*/
        mov r1,r5      /*出力する数をr1へ*/
        ldrpl r0,s1a   /*書式0を使う*/
        ldrmi r0,s0a   /*書式1を使う*/
        bl printf      /*printfで出力*/
        add r5,r5,#1
 cmp r5,#6      /*終りの空白が済んだときのr5の値*/
 bmi l0
 ldr r0,s2a
 bl printf      /*書式2で改行を出力*/
 mov r0,#1      /*プログラムの返り値*/
        ldr r1,reta
        ldr lr,[r1]    /*帰り番地を復活*/
        bx lr

aa:     .word a
s0a:    .word s0
s1a:    .word s1
s2a:    .word s2
reta:   .word ret
想定通りの出力が得られる. (空白の様子が分かるように.や*が入っている.)
. *. *. *. 1. 2. 3. *. *
上のプログラムのミソは
l0:     cmp r5,#1      /*出力する最初の数a*/
        rsbgts r0,r5,#3/*出力する最後の数b*/
        mov r1,r5      /*出力する数をr1へ*/
        ldrpl r0,s1a   /*書式0を使う*/
        ldrmi r0,s0a   /*書式1を使う*/
        bl printf      /*printfで出力*/
        add r5,r5,#1
 cmp r5,#6
 bmi l0
のところだ. ラベルがl0:の行ではr5と1を比較, つまりr5から1を引く. r5が-2から0までは結果は負である. 次の行はこの結果が>0ならrsbつまり逆減算で, 3からr5を引き, 最後のsで結果をCPSRに置く. r5が3を超えると負になるわけだ.

次の行でr5を出力パラメータとしてr1へ置き, 前の計算の正負に従って書式をr0に置き, printfを呼ぶ.

r5は最後の数3を超えても増え続け, つまり空白を出力し続け, 6になるとcmpの結果が負ではなくなり, ループから抜ける.

要するに初めの空白部は, r5との比較が負で判定し, 終りの空白はr5と3を逆に引くことで, 範囲外をともに負にしている.

こんなことが出来るのはARMの命令が多様だからであった.

2015年8月16日日曜日

数式お絵描き

前回のブログでこうもりの描き方がすっかり分ったかといえば, まだすっきりしない部分がある. 下の線を描くプログラムの11から15行目が思いのほか複雑なのだ.

この式の構造はこうなっている. 11行目の(1/2)*に掛けられる因子は4つの項で出来ている.
3*Sqrt[1-(x/7)^2]            (楕円)
+Sqrt[1-(Abs[Abs[x]-2]-1)^2] (4つの小楕円)
+Abs[x/2]                    (折線)
-((3*Sqrt[33]-7)/112)*x^2    (放物線)
-3
最初はお馴染の楕円で高さは3にしてある. 次が4個の楕円で高さはどれも1だ. その次が(0,0)で折り返す直線で, 最後がx2の抛物線である.



4個の楕円はx=-4,-2,0,2,4の5個所でy=0である. このうち中央のx=0では, 他の項の楕円の3を足し, 最後の3を引き, 折線と放物線はx=0で0だから, 0のままである.

x=4では楕円は3√33/7, 折線は2, 放物線は16*(3√33-7)/112, それから3を引くが, これは計算すると結局0になる.

x=2のところは面倒になったので, Schemeで計算した.
(define (a x)
 (define (sq x) (* x x))
 (+ (* 3 (sqrt (- 1 (sq (/ x 7)))))
    (abs (/ x 2))
    (- (* (/ (- (* 3 (sqrt 33)) 7) 112) (sq x)))
    (- 3)))

(a 2) => .5094556875135123
(a 3) => .3881737849967646
(a 1) => .3778577420858067
という次第で, 0,±4では0, ±2では0.5くらいになる. 1と3では1.37程度になろうか.



この図を見ると高さの値は微妙に違うようだが, まぁいいか. この図では±4の外側では定義されていないということで, 青線がないが, 前回のフィルタを掛けた下の図では

のように外側の方もよくみるとx軸上に青の線が見える. 高さも2倍になっているのが分かる. これから楕円の下半分を引くと, 思ったとおりの下の曲線が出来る.

ところでこの下の図は, 前回のこうもりの絵とはすこし違っている. 最後のAspectRatio->Automaticをつけなかったのである.



AspectRatioの詳細はまだわからないが, Plot関数では, AspectRatio->Automatic がないと縦横比が黄金比のφになるらしい.

ディスプレイ上で測ってみると, 左右は56ミリ, 上下は33ミリであった.
(/ 56. 33) => 1.696969696969697
確かに黄金比に近いといえなくはない. x軸の0から2と, y軸の0から2の長さは明らかに違う. Automaticを指定した前回のこうもりの絵では, 横:縦は正しく7:3になっている.

黄金比にするような見栄えを工夫するのがAutomaticの機能のように思うのだが, Mathematicaの文化はこういうものか. 前回のSinの図でも, Automaticがないので, 横/縦は1.7くらい. Automaticを追加すると, 縦と横の長さは1:1になり, もっと平べったい, 正確になった. (世の中のSinカーブでは, 前回の図も含めてx=0での勾配が1になっていない図がかなり多い.)

2015年8月5日水曜日

数式お絵描き

Raspberry PiでMathematicaが使えるのはうれしい.

そのアナウンスにこうもりのような絵があった. もちろんMathematicaのPlot関数で描いたものだ.



その描画プログラムも添えてあったので, Mathematica屋の絵を描く方法が分るのではないかと期待してプログラムを解読した. とりあえずは目がくらくらするが,

である.

そもそもMathematicaで標準関数の絵を描くのは簡単だ. "Plot["と書き, 描くべき式を書き, ","で区切り, "{x,-Pi,Pi}"のように変数の範囲を書き, "]"で閉じる.
Plot[Sin[x],{x,-Pi,Pi}]


こうもりの絵のようなのは, 上の線と下の線を{,}で囲み区切ればよい.

Plot[{Sin[x],Cos[x]},{x,-Pi,Pi}]


こうもりプログラムでは, 1行目の"With"から10行目の"]]"までが上の曲線, 11行目から15行目の"2]"までが下の曲線である.

上の曲線のプログラムは"With[{"から w(3行目), l(2,3行目), h(4,5,6行目), r(7,8行目)を局所定義し(8行目の"},"で終る), それらを使い9,10行目の本体で描画する.

w(3*Sqrt[1-(x/7)^2])はこういう曲線だ.

これは要するに長軸(半径)7, 短軸1のx2/49+y2=1なる楕円の, 平方根は正をとるからその上半分を描いたわけである. あとから縦軸を3倍した.

次はl.(6/7)*Sqrt[10]+(3+x)/2-(3/7)*Sqrt[10]*Sqrt[4-(x+1)^2] のような形だ. 下に左右反対のrがある. xで変る項は(3+x)2とSqrt[4-(x+1)^2]で, 前者は斜線, 後者は例によって楕円である.

大きさを合せて合成するとこうなる. 上の図の斜線はxの範囲-7から7まであるが, 下の図では引くべき楕円が定義されている範囲が-3から1までなので, その外の 線は消えている.



4行目のhはこうもりの耳のところで, ( (1/2)*(3*(Abs[x-1/2]+Abs[x+1/2]+6)-11*(Abs[x-3/4]+Abs[x+3/4])))といろいろの斜線を組合せているに過ぎない.

ここまで準備したところで, 上の曲線は-7から-3,-3から-1,-1から1,1から3,3から7の区間をUnitStepで接続する. UnitStep[x]はx<0では0, 0≤xでは1になる関数で, UnitStep[x+3]とすると, -3から右が1になる.

w+(l-w)*UnitStep[x+3]+(h-l)*UnitStep[x+1]+(r-h)*UnitStep[x-1]+ (w-r)*UnitStep[x-3]]

x=-7では4個のどのUnitStepも0だから最初の項wで描きだす. x=-3になるとUnitStep[x+3]が1になり, それに(l-w)が掛るから, wが消えlの曲線を描き始める. x=-1 でlが終りhが始まる. x=1でrに, x=3で(w-r)があるので, rも終り, wが復活して最後を描く.


下の線もなかなか凝っているが, 全体的には上と同じ楕円に, 4個の楕円が組み合わさっている構造だ.
ここで一番主要な線はSqrt[1-(Abs[Abs[x]-2]-1)^2]だ.



すぐ上の図は左からAbs[x], Abs[x]-2, Abs[Abs[x]-2], Abs[Abs[x]-2]-1で, ここでそれぞれの斜線部分に対応して例の楕円が4個できる仕掛けである.



UnitStepの代わりにこちらではバンド幅フィルタのような関数でこの4個の楕円を取り出している. それが((x+4)/Abs[x+4]-(x-4)/Abs[x-4])で, 左がx<-4で-1, -4<xで1, 右がx<4で-1, 4<xで1だから, 結局x<-4で0, -4<x<4で2, 4<xで0という形が得られる.



これで4個の楕円を切り出しておき, それから最初と同じ楕円を引いて下の線が描ける.



勿論面倒なのでは区間の切れ目で位置を合わせることで, Sqrt[33]などはそういう計算の結果であろうが, プログラムの骨組みの説明では無視させて貰った. 面白かった.

2015年7月16日木曜日

SATソルバ

Algorithm7222Dも動くようになった. 充足可能なのには可能にする値を, 不可能なのには`unsatisfiable'を出力する.

こんなにgoto文だらけのプログラムがよく動くと感心するが, それがこの本のアルゴリズムの特徴だ. とにかくプログラムが動くようになれば, その機能が理解できるという期待から, 私はいつも実装に心掛けている.

今回お目にかけるこのプログラムも, 中に説明は書いてないが, TAOCPの記述に忠実に実装してあるので, 本文と読み較べてほしい. もっともその方の説明も分り難いが.

プログラムリストについているのは, Waerden(3,3;9)用のデータである.

(define (and1 n)
 (if (even? n) 0 1))
(define (xor1 n)
 ((if (even? n) + -) n 1))
(define (rsh1 n)
 (fix:lsh n -1))
(define (iverson b) (if b 1 0))

;Waerden(3,3;9)用のデータ
(define n 9) (define m 32) ;waerden(3,3;9)
(define ls '(3 11 19 7 13 19 5 11 17 3 9 15 11 15 19 9 13 17
7 11 15 5 9 13 3 7 11 15 17 19 13 15 17 11 13 15 9 11 13 7 9
11 5 7 9 3 5 7 2 10 18 6 12 18 4 10 16 2 8 14 10 14 18 8 12 
16 6 10 14 4 8 12 2 6 10 14 16 18 12 14 16 10 12 14 8 10 12 
6 8 10 4 6 8 2 4 6) )
(define ws '
(0 0 1 17 2 18 3 19 4 20 5 21 6 22 7 23 0 0 0 0))
(define link '
(0 8 9 10 11 12 0 0 13 14 15 0 0 16 0 0 0 24 25 26 27 28 0 0
29 30 31 0 0 32 0 0 0))
(define start '(96 93 90 87 84 81 78 75 72 69 66 63 60 57 54
51 48 45 42 39 36 33 30 27 24 21 18 15 12 9 6 3 0))

;途中の様子を出力する関数
(define (printnext next h t)
 (let ((ring (list)))
  (define (ploop)
   (if (eq? h t) (display (reverse (cons t ring)))
    (begin (set! ring (cons h ring))
     (set! h (list-ref next h)) (ploop))))
  (ploop)))
(define (printcl j)
 (define (pin x)
  (string->symbol (string-append
   (number->string (quotient x 2))
   (list-ref '("+" "-") (and1 x)))))
 (display (map (lambda (x) (pin (list-ref ls x))) 
  (a2b (list-ref start j) (list-ref start (- j 1))))))
(define (printx xs) 
 (display (map (lambda (x) (if (< x 0) '- x)) xs)))

(define d 0) (define h 0) (define t 0) (define f 0) 
(define k 0) (define p 0) (define b 0) (define l 0) 
(define j 0) (define j1 0) (define i 0) (define l1 0) 
(define q 0)
(define xs (map (lambda (x) 0) (a2b 0 (+ n 1))))
(define next (map (lambda (x) 0) (a2b 0 (+ n 1))))
(define hs (map (lambda (x) 0) (a2b 0 (+ n 1))))
(define ms (map (lambda (x) 0) (a2b 0 (+ m 1))))

(call-with-current-continuation
 (lambda (exit)

(define (d1)
 (list-set! ms 0 0)
 (set! d 0)
 (set! h 0)
 (set! t 0)
 (do ((k n (- k 1))) ((= k 0)) 
  (list-set! xs k -1)
  (if (or (not (= (list-ref ws (* 2 k)) 0))
          (not (= (list-ref ws (+ (* 2 k) 1)) 0)))
  (begin
   (list-set! next k h) (set! h k)
   (if (= t 0) (set! t k)))))
 (if (not (= t 0)) (list-set! next t h))
(d2))

(define (d2) 
 (newline) (if (not (= t 0)) (printnext next h t)) 
 (printx (cdr xs)) ;xsを出力
 (if (= t 0) (exit 'satisfiable)
  (begin (set! k t) (d3))))

(define (ex135 l)
 (let ((j (list-ref ws l)))
 (define (un0)
  (set! p (+ (list-ref start j) 1))
   (un1))
 (define (un1)
  (if (= p (list-ref start (- j 1))) 1 (un2)))
 (define (un2)
  (if (= (list-ref xs (rsh1 (list-ref ls p)))
         (and1 (list-ref ls p)))
   (begin (set! p (+ p 1)) (un1)) (un3)))
 (define (un3)
  (set! j (list-ref link j))
  (if (= j 0) 0 (un0)))
 (if (not (= j 0)) (un0) 0)))

(define (d3)
 (set! h (list-ref next k))
 (let ((f0 (ex135 (* 2 h))) (f1 (ex135 (+ (* 2 h) 1))))
 (set! f (+  f0 (* 2 f1))) )
 (cond ((= f 3) (display (list 'un 
   (string-append (number->string h) "+")
   (string-append (number->string h) "-"))) (d7))
       ((or (= f 1) (= f 2))
        (display (list 'un (string-append (number->string h) 
         (list-ref '(0 "+" "-") f)))) 
        (list-set! ms (+ d 1) (+ f 3)) (set! t k) (d5))
       ((not (= h t)) (set! k h) (d3))
       (else (d4))))

(define (d4)
 (set! h (list-ref next t))
 (list-set! ms (+ d 1)
  (iverson (or (= (list-ref ws (* 2 h)) 0)
               (not (= (list-ref ws (+ (* 2 h) 1)) 0)))))
 (d5))

(define (d5)
 (set! d (+ d 1))
 (list-set! hs d h)
 (set! k h)
 (if (= t k) (set! t 0)
   (begin (list-set! next t (list-ref next k))
    (set! h (list-ref next k))))
 (d6))

(define (ex136)
  (define (ud0)
   (set! j1 (list-ref link j))
   (set! i (list-ref start j))
   (set! p (+ i 1))
   (ud1))
  (define (ud1)
   (if (= (list-ref xs (rsh1 (list-ref ls p)))
          (and1 (list-ref ls p)))
    (begin (set! p (+ p 1)) (ud1))
    (ud2)))
  (define (ud2)
   (set! l1 (list-ref ls p))
   (list-set! ls p l)
   (list-set! ls i l1)
   (ud3))  
  (define (ud3)
   (set! p (list-ref ws l1))
   (set! q (list-ref ws (xor1 l1))) 
   (if (or (not (= p 0)) (not (= q 0)) 
       (>= (list-ref xs (rsh1 l1)) 0))
    (ud5) (ud4)))
  (define (ud4)
   (if (= t 0) (begin (set! t (rsh1 l1)) (set! h (rsh1 l1))
      (list-set! next t h))
    (begin (list-set! next (rsh1 l1) h) (set! h (rsh1 l1))
      (list-set! next t h)))
    (ud5))
  (define (ud5)
   (list-set! link j p) 
   (list-set! ws l1 j)
   (ud6))
  (define (ud6)
   (printcl j)
   (set! j j1)
   (if (not (= j 0)) (ud0)))
 (set! l (+ (* 2 k) b))
 (set! j (list-ref ws l))
 (list-set! ws l 0)
 (if (not (= j 0)) (ud0)))

(define (d6)
 (set! b (modulo (+ (list-ref ms d) 1) 2))
 (list-set! xs k b)
 (ex136) 
 (d2))

(define (d7loop)
 (if (>= (list-ref ms d) 2) (begin 
  (set! k (list-ref hs d))
  (list-set! xs k -1)
  (if (or (not (= (list-ref ws (* 2 k)) 0))
          (not (= (list-ref ws (+ (* 2 k) 1)) 0)))
   (begin (list-set! next k h)
    (set! h k)
    (list-set! next t h)))
  (set! d (- d 1)) (d7loop))))
   
(define (d7)
 (set! t k)
 (d7loop)
 (d8))

(define (d8)
 (if (> d 0) (begin
  (list-set! ms d (- 3 (list-ref ms d)))
  (set! k (list-ref hs d)) (d6))
 (exit 'unsatisfiable)))

(d1)))
このアルゴリズムの出力を, TAOCPにあるのと同じ場所まで書いてみるとこうなる. (プログラムの実際の出力はもっと冗長だが, 不要な空白など手でけしてある.) 左からアクティブリング, xの値, unはユニットになったリテラル(これが同じ変数の正負のリテラルだとバックトラックに入る,) 書き換えた 項を示す.
(1234567) - - - - - - - - - 2+1+3+ 3+1+5+ 4+1+7+ 5+1+9+
(234567)  0 - - - - - - - - 3+1+2+ 3+2+4+ 4+2+6+ 5+2+8+
(34567)   0 0 - - - - - - -(un 3+) 4-3-5- 5-3-7- 6-3-9-
(4567)    0 0 1 - - - - - - 6+2+4+ 7+1+4+ 5+4+6+ 6+4+8+
(567)     0 0 1 0 - - - - -(un 6+) 9-3-6- 7-6-8-
(975)     0 0 1 0 - 1 - - -(un 9-)
(75)      0 0 1 0 - 1 - - 0(un 7+) 8-6-7- 8-7-9-
(85)      0 0 1 0 - 1 1 - 0(un 8-)
(5)       0 0 1 0 - 1 1 0 0(un 5+5-) 5-3-4- 5-4-6- 6-4-8-
(69785)   0 0 1 1 - - - - -(un 5-) 4+5+6+ 8+2+5+ 9+1+5+ 
                                   6+5+7+ 7+5+9+
(6978)    0 0 1 1 0 - - - -(un 9+) 6-3-9-
3行目の2番目の項はTAOCPのと違う. もちろん向こうが違っている. クヌース先生には連絡済みだ.

2015年7月2日木曜日

SATソルバ

アルゴリズムBでLangford(3)を実行したのが次だ.

Langford(3)だから, 1,2,3を2個ずつ, 例えば323121のように並べる. 但し2個のnの間には他の数がn個なければならないという制限があるから, 先ほどのは駄目で, 231213がひとつの解である.

このくらいならちょっとやっただけで出来るが, 長くなる(nが増える)と大事だ.

TAOCPでの方法はこうだ. まずこういう表を作る.


1行目 1を1番と3番に置く. 間隔は1だ. 2行目 1を2番と4番に置く. ...

5行目 2を1番と4番に置く. ...

8行目 3を1番と5番に置く.

3を2番と6番に置くというのもありだが, 対称だから8行目だけにする.

表が出来たらd1の列に1のあるxjから1つ選ぶ. つまりx1, x2, x3, x4から1つ選ぶ. d2からもx5, x6, x7から1つ選ぶ. これをs6の列まで行う.

1つ選ぶ式をTAOCPではS1と表記するから, この関係は

の式になり, S1の展開式を当てはめると, SAT


が得られる. 2-4-, 1-3- が重複しているからそれを除き, リテラルの内部表現にすると
(define sat '((2 4 6 8) (3 5) (3 7) (3 9) (5 7) (5 9) (7 9)
(10 12 14) (11 13) (11 15) (13 15) (16) 
(2 10 16) (3 11) (3 17) (11 17) (4 12) (5 13) (2 6 14) (3 15)
 (7 15) 
(4 8 10) (5 11) (9 11) (6 12 16) (7 13) (7 17) (13 17) (8 14)
 (9 15)))
となる. 従って
(define ls '(9 15 8 14 13 17 7 17 7 13 6 12 16 9 11 5 11 4 8 
10 7 15 3 15 2 6 14 5 13 4 12 11 17 3 17 3 11 2 10 16 16 13
 15 11 15 11 13 10 12 14 7 9 5 9 5 7 3 9 3 7 3 5 2 4 6 8))
前回と同様にsatの先頭だけとると, sが
(0 2 3 3 3 5 5 7 10 11 11 13 16 2 3 3 11 4 5 2 3 7 4 5 9 
 6 7 7 13 8 9)
になり, zも作れ, wやlinkも得られる.
((0) () (19 13 1) (20 15 14 4 3 2) (22 17) (23 18 6 5) 
(25) (27 26 21 7) (29) (30 24) (8) (16 10 9) () (28 11) 
() () (12) ())

ws=
(0 0 1 2 17 5 25 7 29 24 8 9 0 11 0 0 12 0)

link=
(0 13 3 4 14 6 18 21 0 10 16 28 0 19 15 20 0 22 23 0 0 26 0 0
 30 0 27 0 0 0 0)
節の長さがばらばらだから, STARTも作らなければならない.
(define a (map length sat))
a => (4 2 2 2 2 2 2 3 2 2 2 1 3 2 2 2 2 2 3 2 2 3 2 2 3 2
 2 2 2 2)

(define (ac ls)
 (if (null? ls) '(0)
  (let ((a (ac (cdr ls))))
    (cons (+ (car ls) (car a)) a))))

(define start (ac a))
start =>
(66 62 60 58 56 54 52 50 47 45 43 41 40 37 35 33 31 29 27 24
 22 20 17 15 13 10 8 6 4 2 0))
これで準備が完了し, アルゴリズムBを実行すると01000011が得られ, x2,x7,x8を取るのが解と分かり, 泰山鳴動の一幕であった.

2015年7月1日水曜日

SATソルバ

TAOCP 7.2.2.2にあるSATソルバのうち, 今回はAlgorithm 7222Bの話である.

これも結構難しかった. 演習問題128の解答の(7)の式の初期値が頼りであった.

Algorithm Aに較べると, F, B, C, SIZEの配列はなくなり, 代わりにWとLINKが新たに出来た.

SATの問題では, リテラルは節のいろいろな位置に現れるが, このアルゴリズムでは, 節の先頭に現れるものだけについて初期値に設定する. Wはそれぞれのリテラルについて, 最初に現れた節の番号を示す. どの先頭にも現れなければ, 終という意味で0を入れる.

同じリテラルが複数の節の先頭に現れた時は, 続きをLINKで示す. 具体的には下の例を見てほしい.

7.2.2.2の初めの方にSimple Exampleというのがある. 8ビットのx1...x8で3個の0の間隔が等しくなく, 3個の1の間隔も等しくないものを探せというのである.

すぐ下に9ビットの場合のSATの式があるから, それの8ビット版を書けばよいわけで, 実装に適した私流の記法では

1+2+3+,2+3+4+,3+4+5+,4+5+6+,5+6+7+,6+7+8+,
1+3+5+,2+4+6+,3+5+7+,4+6+8+,
1+4+7+,2+5+8+,
1-2-3-,2-3-4-,3-4-5-,4-5-6-,5-6-7-,6-7-8-,
1-3-5-,2-4-6-,3-5-7-,4-6-8-,
1-4-7-,2-5-8-,
つまり8変数, 24節になる. (n=8, m=24)

これをリテラルの内部表現にすると
(define n 8)

(define sat '
((2 4 6) (4 6 8) (6 8 10) (8 10 12) (10 12 14) (12 14 16) 
(2 6 10) (4 8 12) (6 10 14) (8 12 16) (2 8 14) (4 10 16) 
(3 5 7) (5 7 9) (7 9 11) (9 11 13) (11 13 15) (13 15 17) 
(3 7 11) (5 9 13) (7 11 15) (9 13 17) (3 9 15) (5 11 17)))

(define m (length sat)) => 24
配列Lはこれをreverseし, appendしたものである.

先頭の節を1番として, 2から17までのリテラルがどの節の先頭に現れるかをみると (listは0から始まるのでダミーの0をconsしてある)

(define s (cons 0 (map car sat)))
=>
(0 2 4 6 8 10 12 2 4 6 8 2 4 3 5 7 9 11 13 3 5 7 9 3 5)
リテラル2は1番の節に, 4は3番の節に現れる. リテラル l が現れる節番号のリストを l 番目にもつリストzを作る.
(define z (map (lambda (x) '()) (a2b 0 (+ (* n 2) 2))))

z => (() () () () () () () () () () () () () () () () () ())

(for-each (lambda (x y)
 (list-set! z x (cons y (list-ref z x)))) s (a2b 0 (+ m 1)))

z=>
((0) () (11 7 1) (23 19 13) (12 8 2) (24 20 14) (9 3) (21 15)
 (10 4) (22 16) (5) (17) (6) (18) () () () ())
リテラル2は1,7,11の節に, リテラル3は13,19,23の節に現れるということだ. これからWとLINKを作る.
(define ws (map (lambda (x) 0) (a2b 0 (+ (* 2 n) 2))))

(for-each (lambda (x y) 
(if (not (null? x)) (list-set! ws y (car (reverse x)))))
 z (a2b 0 (+ (* 2 n) 2)))

ws => (0 0 1 13 2 14 3 15 4 16 5 17 6 18 0 0 0 0)

(define link (map (lambda (x) 0) (a2b 0 (+ m 1))))

(for-each (lambda (x) (let ((t 0))
(map (lambda (y) (list-set! link y t) (set! t y)) x))) z)

link =>
(0 7 8 9 10 0 0 11 12 0 0 0 0 19 20 21 22 0 0 23 24 0 0 0 0)
この解は00110011, 01011010,01100110,10011001,10100101,11001100だが, このアルゴリズムは最初の解(00110011)しか返さない.
(define (algorithm7222b)
(define (and1 n) (if (even? n) 0 1))
;1とのandをとる
(define (xor1 n) ((if (even? n) + -) n 1))
;1とのxorをとる. リテラルの正負を反転する
(define (rsh1 n) (fix:lsh n -1))
;1ビット右シフトする
(define (? b) (if b 1 0))
;iverson notation [b]=b?1:0

;simple example
(define n 8) (define m 24)
(define sat '
((2 4 6) (4 6 8) (6 8 10) (8 10 12) (10 12 14) (12 14 16) 
(2 6 10) (4 8 12) (6 10 14) (8 12 16) (2 8 14) (4 10 16) 
(3 5 7) (5 7 9) (7 9 11) (9 11 13) (11 13 15) (13 15 17) 
(3 7 11) (5 9 13) (7 11 15) (9 13 17) (3 9 15) (5 11 17)))
(define ls (apply append (reverse sat)))
(define ws '(0 0 1 13 2 14 3 15 4 16 5 17 6 18 0 0 0 0))
(define link '(0 7 8 9 10 0 0 11 12 0 0 0 0 19 20 21 22 0 0 
 23 24 0 0 0 0))
(define start '(72 69 66 63 60 57 54 51 48 45 42 39 36 33 30
 27 24 21 18 15 12 9 6 3 0))
(define ms '(0 0 0 0 0 0 0 0 0))

(define d 0) (define j 0) (define i 0) (define i1 0)
(define j1 0) (define k 0) (define l 0) (define l1 0)

(call-with-current-continuation
 (lambda (exit)

(define (b1) (set! d 1) (b2))
(define (b2) 
 (if (> d n) (begin (do ((d 1 (+ d 1))) ((> d n))
  (display (xor1 (and1 (list-ref ms d))))) (exit 'ok))
 (begin (list-set! ms d (? (or (= (list-ref ws (* 2 d)) 0)
  (not (= (list-ref ws (+ (* 2 d) 1)) 0))))) 
 (set! l (+ (* 2 d) (list-ref ms d))) 
 (b3))))

(define (b3kloop)
 (cond ((= k i1) (list-set! ws (xor1 l) j) (b5))
  ((< k i1) 
  (set! l1 (list-ref ls k))
  (if (or (> (rsh1 l1) d)
    (even? (+ l1 (list-ref ms (rsh1 l1)))))
   (begin
    (list-set! ls i l1)
    (list-set! ls k (xor1 l))
    (list-set! link j (list-ref ws l1))
    (list-set! ws l1 j)
    (set! j j1))
  (begin (set! k (+ k 1)) (b3kloop))))))

(define (b3jloop)
 (if (not (= j 0)) (begin 
  (set! i (list-ref start j))
  (set! i1 (list-ref start (- j 1)))
  (set! j1 (list-ref link j))
  (set! k (+ i 1))
  (b3kloop)
  (b3jloop))))

(define (b3)
 (set! j (list-ref ws (xor1 l)))
 (b3jloop)
 (b4))

(define (b4)
 (list-set! ws (xor1 l) 0)
 (set! d (+ d 1))
 (b2))

(define (b5)
 (if (< (list-ref ms d) 2) (begin
  (list-set! ms d (- 3 (list-ref ms d)))
  (set! l (+ (* 2 d) (and1 (list-ref ms d)))) 
  (b3))
 (b6)))

(define (b6)
 (if (= d 1) (exit "unsatisfiable")
  (begin (set! d (- d 1)) (b5))))

(b1)))

(algorithm7222b) => 00110011
他の実行例については次回にでも.

2015年6月30日火曜日

軌跡作図器

JavaScriptによるアニメーションができた.



ここにあります.

PがO'を中心とする円上を移動すると, Qは直線上を移動する.

2015年6月27日土曜日

軌跡作図器

先日のブログの最後にJavaScriptで動かしたいと書いた. まだそこまではいっていないが, Pの位置を変えるとQがどう動くが調べたくなって, PostScriptでプログラムを書いた.

方針としてはこうだ. AB, AO, ADの長さを決める. O, O'の位置を決める.

OO'の長さに等しくO'Pをとるわけだが, つまりO'を中心とする半径OO'の円上をPが動くわけだが, Oの反対側をO''とし, ∠O''O'Pを0から順に増やすことで, Pの位置を決めることにした.

OとPの位置やOBとPBの長さが分っているので, Bの位置が決まる. Bが決まるとAが決まり, BとPが分るからCも決まる. AとCからの長さが分っているから, Dの位置が決まり, Qも得られる.

このようにして, ∠O''O'Pを0°, 30°, 60°, 90°, 120°, 135°と変えて描いたのが次の図である.





30°


60°


90°


120°


135°

この先Pはどこまで進めるかというと, OBとPBの長さが決まっているから, OPの間隔がPB-OBになったときに一直線になって以後破綻する.

2015年6月26日金曜日

軌跡作図器

學士會会報No.912 (2015-III)の表紙の図に興味をもった. 説明には独国マルチン・シリング(ハーレ社)製「軌跡作図器(Serie XXIV, Modell nr.11)」とあるが, それに続く説明は要領を得ないので改めて説明したい.

同じような図をネットで探して頂いてきたのが次の図だ.(http://goo.gl/Fpfcqt)



リンクABとCDの長さは等しく, ADとBCの長さも等しい. A,B,C,Dはピボットで連結されていて角度は変わる. AB, CB, ADのそれぞれをm : (1-m)に内分した点をO,P,Qとする(0 < m < 1).

Oは固定点で, ABはOを支点に回転する. PにPO'がOO'に等しくなるようなリンクPO'を取付け, O'も固定する. つまりPはO'を中心に円を描き, この円はOを通る.

こうしてPが円(の一部)を描くと, Qは直線の軌跡を描くという道具であった(Hartのinversorという名前である).

上の図を次のように書き直す. 太線がリンク機構である. 三角形ACBと三角形CADは逆向きだが合同だから, ACとBDは平行. ABDCは等脚台形になる. OPQは一直線でこれもACに平行.



BD上にA, CからACに垂線を立ててE,Fとする. EB=DFだから
AC×BDはEF×BD=(ED+EB)×(ED-EB)=ED2-EB2.
ED2+AE2=AD2
EB2+AE2=AB2
だから
AC×BC=AD2-AB2
AO/AB=CP/CB=AQ/AD=mだから
OP=m×AC, OQ=(1-m)×BD
従って
OP×OQ=m×(1-m)×(AD2-AB2)=一定.

OP'×OQ'=OP×OQ=一定 でP'がOPの中点O'を中心とする円を動くと, Q'はQでOQに直交する線上を動く.



上の図で∠OP'P=∠R. ∠PQQ'=∠R だから4点PQQ'P'は円周上にあり, 方冪の定理(定点を過ぎる直線が定円と交わるとき, 定点より2つの交点までの長さの積は一定) によりOP'×OQ'=OP×OQ.

そのうちJavaScriptで動かしたいと思う.

2015年6月22日月曜日

SATソルバ

TAOCP 7.2.2.2にはSATソルバのアルゴリズムがいくつか登場する.

今回はそのうちAlgorithm 7222Aを話題としたい. ただこのアルゴリズムの背景を初めから述べるのは大変なので, TAOCPのその辺を読んである読者を対象に書くことで勘弁して貰いたい.

TAOCPのアルゴリズムは自然言語であるだけでなく, 非常に簡潔な記述なので, 前後の説明をよく理解していないと, なんのことだか分からないことが多い. そういう場合にSchemeで記述したこのブログのプログラムは理解の助けになのではないか. とにかく実際に動くプログラムとして存在するから だ.

TAOCPにも, あるデータを例としての説明はある. 7.2.2.2-6と7の式



で, (6)は充足不能, (7)は充足可能な例である.

まず(7)について, 次のようなデータを用意する.



配列L, F, B, Cは添字が0から30で, Lには(7)の節が置いてある. L(10)から始まる9,7,3は(7)の最後の節の3-,4-,1-を4,3,1の順にし, 正のリテラルは変数番号を2倍し, 負のリテラルは2倍して1足したものなので, 9,7,3になっている. Cの対応する場所の7は, 7番目の節を示す.

L(13)以降も同様である.

FとBは両方向リンクで, F(2)つまりリテラル1+が現れる場所はL(30)であり, F(30)=24は次がL(24)であることを示す. Bは逆向きの情報である.

リンク構造なのは, リテラルや節を削除/挿入するのに便利なためである.

C(2)からC(9)まではそれぞれのリテラルが何個あるかを示す.

それらの下のSTARTは, それぞれの節の先頭の位置; SIZEは節にあるリテラルの個数である.

Schemeで実装したプログラムは以下のようだ.


(define (algorithm7222a)
(define (and1 n) (if (even? n) 0 1))
;1とのandをとる
(define (xor1 n) ((if (even? n) + -) n 1))
;1とのxorをとる. リテラルの正負を反転する
(define (? b) (if b 1 0))
;iverson notation [b]=b?1:0

(define n 4) (define m 7)
(define ls '(0 0 0 0 0 0 0 0 0 0 9 7 3 8 7 5 6 5 3 8 4 3
 8 6 2 9 6 4 7 4 2))
(define fs '(0 0 30 21 29 17 26 28 22 25 9 7 3 8 11 5 6 15 12
 13 4 18 19 16 2 10 23 20 14 27 24))
(define bs '(0 0 24 12 20 15 16 11 13 10 25 14 18 19 28 17 23
 5 21 22 27 3 8 26 30 9 6 29 7 4 2))
(define cs '(0 0 2 3 3 2 3 3 3 2 7 7 7 6 6 6 5 5 5 4 4 4
 3 3 3 2 2 2 1 1 1))
(define size '(0 3 3 3 3 3 3 3))
(define start '(0 28 25 22 19 16 13 10))

(define ms (map (lambda (x) 0) (a2b 0 (+ n 1))))

(define l 0) (define a 0) (define d 0) (define p 0) 
(define i 0) (define j 0) (define q 0) (define r 0) 
(define s 0)

(call-with-current-continuation
 (lambda (exit)

(define (a1)
 (set! a m) (set! d 1)
 (a2))
(define (a2)
 (set! l (* 2 d))
 (if (<= (list-ref cs l) (list-ref cs (+ l 1)))
  (set! l (+ l 1)))
 (list-set! ms d (+ (and1 l) 
  (* 4 (? (= (list-ref cs (xor1 l)) 0)))))
 (newline) (display (list 'ms (take (+ d 1) ms)))
 (if (= (list-ref cs l) a) (begin (newline)
  (do ((i 1 (+ i 1))) ((> i n)) 
   (display (and1 (xor1 (list-ref ms i)))))
  (exit 'satisfiable)))
 (a3))

(define (a3back)
 (if (>= p (+ n n 2)) (begin
  (set! j (list-ref cs p))
  (set! i (list-ref size j))
  (list-set! size j (+ i 1))
  (set! p (list-ref bs p))
  (a3back)) (a5)))

(define (a3loop)
 (if (>= p (+ n n 2)) (begin
 (set! j (list-ref cs p))
 (set! i (list-ref size j))
 (if (> i 1) (begin
   (list-set! size j (- i 1)) 
   (set! p (list-ref fs p)) (a3loop))
  (begin (set! p (list-ref bs p)) (a3back))))
(a4)))

(define (a3)
 (set! p (list-ref fs (xor1 l)))
 (a3loop))

(define (a4loop)
 (if (>= p (+ n n 2)) (begin
  (set! j (list-ref cs p))
  (set! i (list-ref start j))
  (set! p (list-ref fs p)) 
  (do ((s i (+ s 1))) ((= s (+ i (list-ref size j) -1)))
   (set! q (list-ref fs s))
   (set! r (list-ref bs s))
   (list-set! bs q r)
   (list-set! fs r q)
   (list-set! cs (list-ref ls s)
    (- (list-ref cs (list-ref ls s)) 1)))
  (a4loop))))

(define (a4)
 (set! p (list-ref fs l))
 (a4loop)
 (set! a (- a (list-ref cs l))) (set! d (+ d 1)) 
 (a2))

(define (a5)
 (if (< (list-ref ms d) 2) (begin 
  (list-set! ms d (- 3 (list-ref ms d)))
  (set! l (+ (* 2 d) (and1 (list-ref ms d)))) (a3))
 (a6)))

(define (a6)
 (if (= d 1) (exit "unsatisfiable")
  (begin (set! d (- d 1))
   (set! l (+ (* 2 d) (and1 (list-ref ms d)))) (a7))))

(define (a7loop)
 (if (>= p (+ n n 2)) (begin
  (set! j (list-ref cs p))
  (set! i (list-ref start j))
  (set! p (list-ref bs p))
  (do ((s i (+ s 1))) ((= s (+ i (list-ref size j) -1)))
   (set! q (list-ref fs s))
   (set! r (list-ref bs s))
   (list-set! bs q s)
   (list-set! fs r s)
   (list-set! cs (list-ref ls s)
    (+ (list-ref cs (list-ref ls s)) 1)))
 (a7loop))))

(define (a7)
 (set! a (+ a (list-ref cs l)))
 (set! p (list-ref bs l))
 (a7loop)
 (a8))

(define (a8loop)
 (if (>= p (+ n n 2)) (begin
  (set! j (list-ref cs p))
  (set! i (list-ref size j))
  (list-set! size j (+ i 1))
  (set! p (list-ref fs p))
  (a8loop))))

(define (a8)
 (set! p (list-ref fs (xor1 l)))
 (a8loop)
 (a5))

(a1)))
)
(algorithm7222a)
途中, ステップA2で動作状況を示すmsを印字しながらの実行結果は
(ms (0 1))
(ms (0 1 0))
(ms (0 1 0 1))
(ms (0 1 0 1 4))
0101
;... done
;Value: satisfiable
で, 解はx1=0, x2=1, x3=0, x4=1 である.

TAOCPには(6)の実行中のmsの値も掲載してあるから
(define n 4) (define m 8)
(define ls '(0 0 0 0 0 0 0 0 0 0 9 5 2 9 7 3 8 7 5 6 5 3 8 4
 3 8 6 2 9 6 4 7 4 2))
(define fs '(0 0 33 24 32 20 29 31 25 28 9 5 2 10 7 3 8 14 11
 6 18 15 16 4 21 22 19 12 13 26 23 17 30 27))
(define bs '(0 0 12 15 23 11 19 14 16 10 13 18 27 28 17 21 22
 31 20 26 5 24 25 30 3 8 29 33 9 6 32 7 4 2))
(define cs '(0 0 3 3 3 3 3 3 3 3 8 8 8 7 7 7 6 6 6 5 5 5 4 4
 4 3 3 3 2 2 2 1 1 1))
(define size '(0 3 3 3 3 3 3 3 3))
(define start '(0 31 28 25 22 19 16 13 10))
として(6)もやってみると
(ms (0 1))
(ms (0 1 1))
(ms (0 1 1 0))
(ms (0 1 1 3 1))
(ms (0 1 2 1))
(ms (0 1 2 1 1))
(ms (0 1 2 2 1))
(ms (0 2 1))
(ms (0 2 1 1))
(ms (0 2 1 1 1))
(ms (0 2 1 2 1))
(ms (0 2 2 1))
(ms (0 2 2 2 1))
;... done
;Value 11: "unsatisfiable"
のように, TAOCPにあるのと同様な出力が得られた. まぁこの実装は間違ってはいないようだ.

2015年6月1日月曜日

SATソルバ

TAOCP 7.2.2.2は充足可能性(Satisfiability)が話題だ. 最初の方に多くの問題が CNF(Conjunctive Normal Form 乗法標準形)で記述できるという話が続く. 私はこの分野にはさほど興味がなく, まったくのど素人で, いつものようにアルゴリズ ムにだけ関心がある. 繰り返すまでもなく, TAOCPのアルゴリズムは自然言語で 記述されているから, その英文解釈が煩わしく, いつもああでもない, こうでもない と思いつつSchemeに変換し, いくつかの例題が無事に解けるのをみてやっと 安心している.

今回もそういう話である.

早速例をTAOCPから拝借すると, 乗法標準形は次のようなものである.


x1, x2, ...などが変数(variable). 変数または その上に線を引いたものがリテラル(literal). リテラルを∨でつないだものが 論理和節(clause), そしてCNFはそれらの節を∧でつないだものである.

かっこや∧や∨をいちいち書くのも面倒だということで, TAOCPではこれを


のように書いた.

さてSATソルバはFが真になるように, これらの変数に真(1), 偽(0)を対応させる もので, この例のサイズなら目の子で探せば解が見つかるが, 実用になる問題では 変数の個数がすごく多く, 高速なアルゴリズムが研究対象になっているらしい. ( 4月23日のブログMcGregorグラフを解く場合は440変数の1406論理和節になるとか.)

TAOCPにはまず基本的なアルゴリズムが出ていたので, 練習としてそれを Schemeで実装してみた.

まず上のCNFの例は, 私の実装ではS式で((1+ 2-) (2+ 3+) (1- 3-) (1- 2- 3+)) と書くことにする. 負のリテラルに-がつくだけでなく, 正のリテラルに+が つくのは対称性からきれいではないか. Schemeでは+1は数だが, 1+はシンボルになる. 従って 1+1-にしたり, その逆にしたりは
(define (neg l)
(let* ((s (symbol->string l)) (w (- (string-length s) 1)))
 (string-set! s w (integer->char 
  (- 88 (char->integer (string-ref s w)))))
 (string->symbol s)))
つまりリテラルlを貰い, その文字列sとその長さ-1のwを 作る. 最後の+-を取り出し, ASCIIの値とし,
(char->integer #\+) => 43
(char->integer #\-) => 45
だからその和の88から引くと+, -が逆転できる. ついでだが, 変数の記号からリテラルを作る(例えば x からx+, x-, 2 から2+, 2-を作る)手続きは次のようだ.
 (define (lit v s)
  (string->symbol (string-append
   ((if (number? v) number->string symbol->string) v)
   (symbol->string s))))
(lit 'x '+) => x+
(lit 'x '-) => x-
(lit 2 '+) => 2+
(lit 2 '-) => 2-
TAOCPにあるSATソルバの基本的なアルゴリズムB(F)(7222-56)は次のようだ.





考えてみると, このアルゴリズムは当然と思える. Fは節がandで繋がってい るから, 空のandは真だ. (and) => #t.

andで繋がっている節に空のものがあると, 空の節は(or) => #f だから Fは偽になる.

そのどちらでもないなら, まだFに残っている変数の正と偽のリテラル を作り, それで「ただし」以下のF|lを用意してBに 再帰的に作用させ, その結果にリテラルを足したのを返す.

「ただし」の内容はこうだ. lを含む節は真なので, andの並びでは 考慮する必要がないから削除する. 一方lの否定を含む節はorに 偽の項は関係ないから項を削除する. いずれにしろlの変数は 消えることになる.

そうだと思って書いたのが次のSchemeのプログラムである.
(define (b f vs)
 (define (lit v s)
  (string->symbol (string-append
   ((if (number? v) number->string symbol->string) v)
   (symbol->string s))))
 (let ((a '()))
  (define (b1 f ls vs)
   (display (list 'b1 f ls vs )) (newline)
   (cond ((member '() f) #f)
         ((= (length (nub f)) 1)
          (set! a (cons (cons (caar f) ls) a)))
         (else
         (let ((l (lit (car vs) '+)))
           (b1 (apb f l) (cons l ls) (cdr vs)))
         (let ((l (lit (car vs) '-)))
           (b1 (apb f l) (cons l ls) (cdr vs))))))
 (b1 f '() vs) a))
引数のvsは変数のリストである.

だから上の例では
(b '((1+ 2-) (2+ 3+) (1- 3-) (1- 2- 3+)) '(1 2 3))
と呼び出す. 内部手続きb1が途中経過を示すので, 出力は

(b1 ((2+ 3+) (3-) (2- 3+)) (1+) (2 3))
(b1 ((3-) (3+)) (2+ 1+) (3))
(b1 (()) (3+ 2+ 1+) ())
(b1 (()) (3- 2+ 1+) ())
(b1 ((3+) (3-)) (2- 1+) (3))
(b1 (()) (3+ 2- 1+) ())
(b1 (()) (3- 2- 1+) ())
(b1 ((2-) (2+ 3+)) (1-) (2 3))
(b1 (()) (2+ 1-) (3))
(b1 ((3+)) (2- 1-) (3))
=> ((3+ 2- 1-))
トレースの最初の行は, 後ろから見て残りの変数は2と3, リテラルは1+. 従って((1+ 2-) (2+ 3+) (1- 3-) (1- 2- 3+))の左端の項は1+があるから削除. 次は1+, 1-もないからそのまま, 次の2項は1-を削除. そして((2+ 3+) (3-) (2- 3+))を得る.

次の行はリテラルが2+. そこで(2+ 3+)は削除. (2- 3+)の 2-を削除. ((3-) (3+))が残る.

次. 変数3に対して3+, 3-をリテラルとしても()が得られ. εなので 充足不能になる. バックトラックして2-をリテラルにしても失敗.

改めて1-をリテラルにしてみたところ, fが空リストになり成功した.

原理的にはこれでよいが, 実用的には工夫があるらしい. そういう改良 アルゴリズムの話がTAOCPでは続く. その説明はまたいずれ.

2015年4月27日月曜日

Christopher StracheyのGPM

順列の生成


n個の要素のすべての順列を生成するのも情報科学標準問題である.

私のSchemeのライブラリには
(define (permutation ls) ;list of permutation of ls
 (define (list-del l n)
   (if (= n 0) (cdr l)
    (cons (car l) (list-del (cdr l) (- n 1)))))
  (if (null? ls) '(())
   (apply append (map (lambda (i)
    (let ((x (list-ref ls i)))
     (map (lambda (p) (cons x p))
       (permutation (list-del ls i))  )))
      (a2b 0 (length ls))))))
;(permutation '(a b c))=>
;((a b c) (a c b) (b a c) (b c a) (c a b) (c b a))
というのが置いてある.

(list-del l n)はlのn番目の要素を取り去ったリストを返す. 本体は順列をとるリストが空なら空リストのリストを返す. そうでないなら最後の行(a2b 0 (length ls))でlsが(a b c)なら(0 1 2)のリストを作り, その各々を(lambda (i)..)のiとして(let..以下をやったものをappendする.

(let以下のxはlist-refだからa,b,cになり, それぞれlist-delした(b c) (a c) (a b)の順列の先頭に付ける. だから(a b c) (a c b) (b a c) (b c a) (c a b) (c b a)をappendすることになり, (a b c)の順列が完成する.

さてgpmには配列もリストもないから, また引数列を活用するプログラムを考えなければならない. 引数の長さに従って別のマクロを用意しなければならない.

そう考えて書いたのが次だ.
$def,p2,<~1~2~3>;
$def,p3,<$p2,~1,~2,~3;
$p2,~1,~3,~2;>;
$def,p4,<$p3,~1~2,~3,~4;
$p3,~1~3,~2,~4;
$p3,~1~4,~2,~3;>;
$def,p5,<$p4,~1~2,~3,~4,~5;
$p4,~1~3,~2,~4,~5;
$p4,~1~4,~2,~3,~5;
$p4,~1~5,~2,~3,~4;>;
これはまずp5を$p5,,a,b,c,d;のように呼ぶ. 第1引数が空, 第2引数以降がa,b,c,d である.

そこでp5の定義を見ると, $p4,空a,b,c,d;$p4,空b,a,c,d;$p4,空c,a,b,d;$p4,空d,a,b,c;とp4を呼び出している. つまりp4の第1引数は上のSchemeのプログラムのx, 以降はlist-delの結果に対応している.

p4の定義を見ると, $p3,~1~2,~3,~4;と呼ぶから, 最初のp4の呼出からは$p3,空ab,c,d;$p3,空ac,b,d;$p3,空ad,b,c;と呼出され, その最初のp3の呼出では$p2,空ab,c,d;$p2,空ab,d,c;と呼出され, p2によってabcd, abdcと出力される.

なんだか全部の順列を自分で書いたみたいな気もするが, これで完成である. 引数の数だけマクロを用意するのは面倒だが, それを我慢すれば存外簡単であった.

2015年4月26日日曜日

Christopher StracheyのGPM

クィーンパズル生成マクロ


前回は4クィーン問題をGPMで解いてみたが, そもそもは8クィーンを解きたい.

すこし時間が経ったのでおさらいするとして, q2とq3を並べて見る.

$def,q3,<$z,0,                  $def,q2,<$y,0,                  
 $def,z,<$~1,                    $def,y,<$~1,                   
 $def,~1,<$                      $def,~1,<$                     
  $|,$?,>>~1<<,>~1<,0;,           $|,$?,>>~1<<,>~1<,0;,         
  $|,$?,>>~1<<,>~1<,3;,           $|,$?,>>~1<<,>~1<,2;,         
  $|,$?,>>~2<<,>~1<,0;,           $|,$?,>>~2<<,>~1<,0;,         
  $|,$?,>>~2<<,>~1<,2;,           $?,>>~2<<,>~1<,1;;;;,         
  $|,$?,>>~3<<,>~1<,0;,                                         
  $?,>>~3<<,>~1<,1;;;;;;,                                       
  $def,f,<$q4,>>>~1<<<,>>>~2<<<,  $def,f,<$q3,>>>~1<<<,>>>~2<<<,
   >>>~3<<<,>>~1<<,>>>~4<<<;       >>~1<<,>>>~3<<<;             
   $z,$1+,>>~1<<;;>;,              $y,$1+,>>~1<<;;>;,           
  $def,t,<$z,$1+,>>~1<<;;>;;>;,   $def,t,<$y,$1+,>>~1<<;;>;;>;, 
 $def,>~4<,;;>;;>;               $def,>~3<,;;>;;>;              
q2は$q2,0,2,4; $q2,0,3,4; のように既に2個のクィーンが無事に置けたとしてそれらのクィーンの位置と4クィーンであるという4をもって呼ばれる. 3個めを0,1,2,3と置いてみて, うまく置ければそれを次の引数に追加してq3を呼ぶ仕掛である.

3個めを置くマクロはyで, $y,0;と呼ぶ. 続けてyの定義がある. まず引数が~3, つまりnクィーンのnなら終りのチェックをするが, 例によってelseから書いてある. else節はor(|)で$?をいくつも呼んで, 引数で持ち込まれたクィーンと当るかを調べる.

$?,>>~1<<,>~1<,0;は最初のクィーンと今度のクィーンが横並びか見る. (>>~1<<が最初のクィーン, >~1<がyが作った新しいクィーン.) 次の?はクィーンが横方向に2ずれているから, 高さ方向の差が2であるか見る.

次は2番目のクィーンと今度のクィーンが横並びか高さの差が1かを見る.

orが偽をいうことは, 無事に置けたのだから, 受け取った引数と, 今回テスト中のクィーンとnの値4をもってq3を呼び, 今回のクィーンの位置をひとつ増してyを呼ぶ.

新しい位置が当っていれば, $def,t,にあるように, 位置を増してyを呼ぶ.

新しいクィーンの位置がnまで来ると最後の方の$def,>~3<,で, そのまま脱出してしまう.

これが別れば一般のクィーンの個数のqマクロが書ける.

しかしマクロを生成するマクロには宿命がある. 評価を阻止する<と>について, 出力マクロに入れるものか, 現状の評価中に評価を阻止するためのものか, 区別が必要なことである.

今回は出力に書き出すものは, [と]を使うことにした. 出来上がったマクロをエディタで処理し, [,]を<,>に置き換える.

$def,g0,<$~1,
$def,~1,<
<  $|,$?,]]~>$-,>~2<,>~1<;<[[,]~1[,0;,>
<  $|,$?,]]~>$-,>~2<,>~1<;<[[,]~1[,>>~1<<;,>
    $g0,$1-,>~1<;,>~2<;>;,
$def,1,<
<  $|,$?,]]~>$-,>~2<,>~1<;<[[,]~1[,0;,>
<  $?,]]~>$-,>~2<,>~1<;<[[,]~1[,>>~1;;>;
$def,g1,<$~1,
$def,~1,<;;$g1,$1-,>~1<;;>;,
$def,0,;;>;
$def,g2,<$~1,
$def,~1,<<]]]~>$-,>~2<,>~1<;<[[[,>
   $g2,$1-,>~1<;,>~2<;>;,
$def,0,;;>;

$def,gen,<
<$def,q>~1<,[$>~2<,0,>
< $def,>~2<,[$~1,>
< $def,~1,[$>$g0,~1,$1+,~1;;$g1,~1;<,>
<   $def,f,[$q>$1+,~1;<,>
  $g2,~1,$1+,~1;;<]]~1[[,]]]~>$1+,~1;<[[[;>
<   $>~2<,$1+,]]~1[[;;];,>
<  $def,t,[
$>~2<,$1+,]]~1[[;;];;];,>
< $def,]~>$1+,~1;<[,;;];;];>
>;
この$genが生成マクロである. $gen,3,z;のように呼ぶ. するとq3が出力される. zは上のq2の例にあったyの働きのものである.

真中より少し下のgenの定義を見ると,

<$def,q>~1<,[$>~2<,0,>
< $def,>~2<,[$~1,>
< $def,~1,[$>
とあって, これで
$def,q3,[$z,0,
 $def,z,[$~1,
 $def,~1,[$
までが出来る. 次のg0はorの連続を出力, g1はその後のセミコロンの列を出力する. さらに先のg2はq4の呼出し列を生成する. そういう次第で出力されたq3は
$def,q3,[$z,0,
 $def,z,[$~1,
 $def,~1,[$
  $|,$?,]]~1[[,]~1[,0;,
  $|,$?,]]~1[[,]~1[,3;,
  $|,$?,]]~2[[,]~1[,0;,
  $|,$?,]]~2[[,]~1[,2;,
  $|,$?,]]~3[[,]~1[,0;,
  $?,]]~3[[,]~1[,1;;;;;;,
   $def,f,[$q4,]]]~1[[[,]]]~2[[[,]]]~3[[[,]]~1[[,]]]~4[[[;
   $z,$1+,]]~1[[;;];,
  $def,t,[$z,$1+,]]~1[[;;];;];,
 $def,]~4[,;;];;];
ほかのq1,q2,...,q7も以下の呼出しで生成される.
$gen,7,d;
$gen,6,c;
$gen,5,b;
$gen,4,a;
$gen,3,z;
$gen,2,y;
$gen,1,x;
その出力のかっこを変換し, 両端の
$def,q8,<~1,~2,~3,~4,~5,~6,~7,~8;>;
$def,q0,<$w,0,
 $def,w,<$~1,
 $def,~1,<$q1,>~1<,>>~1<<;
  $w,$1+,>~1<;;>;,
 $def,>~1<,;;>;;>;
を追加し, $q0,8; で呼ぶと, 約28分の実行の後
0,4,7,5,2,6,1,3;0,5,7,2,6,3,1,4;
0,6,3,5,7,1,4,2;0,6,4,7,1,3,5,2;
1,3,5,7,2,0,6,4;...
...
7,2,0,5,1,4,6,3;7,3,0,2,5,1,6,4;
が得られた. 結構大騒ぎだが, もう何クィーンでも対応できるようになった. もっとも引数の個数の制限は別である.

2015年4月23日木曜日

McGregorグラフ

TAOCPに「10次のMcGregorグラフ(McGregor graph of order 10)」という図があった. これには「Scientific Americanの1975年4月号で, Martin Gardnerがこの図の塗り分けには5色が必要と紹介し世間を驚かせた」と説明がある.



その号のMathematical Gamesコラムに はeπ√163は整数である(これについては高橋秀俊先生が1975年の高橋コンファレンスで話された. どこかにそれが書かれているかと探したが見当たらない.)とか怪しい話もあり, 結局は四月馬鹿の話題であったのだ.

私もその号のコラムを読み, なーんだと思ったひとりである. それから既に40年経った.

TAOCPには演習問題の解を見る前に自分で試みよとあったので, とりあえずバックトラックしながら解を探すプログラムをSchemeで書き, 解をひとつ見つけた. 結構時間が掛ったが, プログラムを走らせてから長めのミーティングに出ていて, 戻ってきたら解が出力されていた.

ちゃんと4色で塗れた. 探索プログラムより区画どうしの接続データを作るのが面倒であった. またこの図を描くPostScriptのプログラムも予想外に時間がかかった.

接続データはこういう形である.
(0 1 10 11 100 101 102 103 104 105 109)
(1 0 2 11 12 109)
(2 1 3 12 13 109)
(3 2 4 13 14 109)
(4 3 5 14 15 109)
(5 4 6 15 16 109)
... 途中100行省略
(106 10 95 96 105 107)
(107 10 96 97 106 108)
(108 10 97 98 107 109)
(109 0 1 2 3 4 5 6 7 8 9 10 9 19 98 108)))
例えば最上行は「0の区画と隣合うのは1, 10, 11, a0, a1, a2, a3, a4, a5, a9である」という意味. aは10 と書いてある. もちろんaの隣のリストにあるbについて, bの隣のリストにaがあることはチェックした.

得られた解で塗り分けたのが下の図である. かなり美しい.



110の区画でそれぞれの色が使われた数はであった. と書いた理由はTAOCPにひとつの色は7区画でしか使わない解があるという記述があったからだ. そういう解を探すのはまた大変そうだ.

TAOCPではこれをdancing linksで解くつもりらしい. dancing linksのプログラムは以前に書いたこともあるので, そのうちdancing linksでもプログラムを書いてみたい.

ところで各区画に付いている番号には意味があるのかな. 上の図は10次だが, 9次とか描くには番号が手引きになるのだろうか.