2010年5月8日土曜日

Life Game

この所 TAOCPのアルゴリズムの解説ブログの気味があるが....

John Conwayが発明したとされるLife Gameは, 1970年にMartin GardnerがScientific American誌で紹介し, 爆発的に広まった. 私もそのScientific American誌のコラムを読んだ時は興奮した. その時からしばらく, 多くの計算時間がグライダーやグライダーガンの発見に費やされたという.

要するにMoore近傍, 状態0と1のセルオートマトンで, 中央のセルが0の時, 周囲の和が3なら1, それ以外は0. 中央が1の時, 周囲の和が2か3なら1, それ以外は0という遷移規則である.もっとも研究されたセルオートマトンである.

TAOCPでのLife Gameの記述に興味がある. まず遷移規則は(演習問題7.1.3-167)
f(xnw,...,x,...,xse)=[2<xnw+xn+xne+xw+1/2x+xe+xsw+xs+xse<4]
である. xは中央のセル, xnwはxの北西のセルのように読む.
[boole式]は, Iverson bracketで, 中のboole式が真なら1, 偽なら0を表す. xが0なら, 和が3の時だけ真だ. xが1なら, 1/2があるので, 和が2.5と3.5の時真になる; つまり周囲の和が2か3ということだ.

この遷移式fを26個のboole連鎖で表せるという. どうするか. 加算するには全加算器が必要. 全加算器
(z1z0)2=x1+x2+x3
は5個のboole演算で実現出来る. すなわち
x4=x1⊕ x2,
z0=x5=x3⊕ x4;
x6=x3∧ x4,
x7=x2∧ x2,
z1=x6∨ x7.
z0はx1, x2, x3の排他論理和だからこうであろう. z1は, xの2つが1なら1なので, x7はx1x2が共に1の時, x4はx1とx2のいずれかが1だから, x6はとx3とx1かx2が共に1の時で, その論理和でz1が得られる.

さて
(z1z2)2=xnw+xw+xsw,
(z3z4)2=xn+xs,
(z5z6)2=xne+xe+xse,
(z7z8)2=z2+z4+z6,
f=S1(z1,z3,z5,z7)∧(x∨z8),
where
S1(x1,x2,x3,x4)=((x1∨x2)∨(x3∧x4))⊕((x1∧x2)∨(x3∨x4)).
でfが計算出来る. boole演算の個数は全加算器(5x3), 半加算器(2), S1(7), それに最後の∧と∨で合計26になる.

これでfが計算出来る理由は次のようだ.
(z1z2)2=xnw+xw+xswからz2は左の列のパリティである. z1は左の列の1が2個か3個あることを示す. 同様にz4は, 中の列の上下のセルのパリティ, z3は, 中の列に1が2個あることを示す. z6とz5についても同様である.

(z7z8)2=z2+z4+z6から, z8はxを除いた全体のパリティである. またz7は, 各列のパリティの和が2か3かを示す.

ところで対称関数S1は, 4個の引数のうち, 1が1個の時に限り1を返す. 従ってS1(z1,z3,z5,z7)が1になるのは,
1. 左の列だけに1が2個か3個あるか,
2. 中の列に1が2個あるか,
3. 右の列に1が2個か3個あるか,
4. 列のパリティの和が2か3かのいずれかの時である.
1の場合なら, 左に1が2個か3個あるだけで, 他は0である.
2なら, 中に1が2個あるだけ, 3なら右の列に2個か3個あるだけだ.
4の場合はどうか. 他の列の1の数は0個か1個である. パリティの和が2以上だから, 1が1個の列が2列か3列あるわけで, やはり1の数は全体で2個か3個である.

xが1なら, 周りに2個か3個の1があればxは生きるからこれでOKだ. xが0なら, z8のパリティが1なら, 周りに3個の1があることになり, xは1になれる.
これがこの計算のからくりであった.

私のSchemeでの実装は次の通り.

(define (f xnw xn xne xw x xe xsw xs xse)
(define (and x y) (* x y))
(define (or x y) (quotient (+ x y 1) 2))
(define (xor x y) (modulo (+ x y) 2))
(define (fa x1 x2 x3)
(let* ((x4 (xor x1 x2)) (x6 (and x3 x4)) (x7 (and x1 x2)))
(cons (or x6 x7) (xor x3 x4))))
(define (ha x1 x2)
(cons (and x1 x2) (xor x1 x2)))
(define (s1 x1 x2 x3 x4)
(xor (or (or x1 x2) (and x3 x4))
(or (and x1 x2) (or x3 x4))))
(let* ((z12 (fa xnw xw xsw)) (z1 (car z12)) (z2 (cdr z12))
(z34 (ha xn xs)) (z3 (car z34)) (z4 (cdr z34))
(z56 (fa xne xe xse)) (z5 (car z56)) (z6 (cdr z56))
(z78 (fa z2 z4 z6)) (z7 (car z78)) (z8 (cdr z78)))
(and (s1 z1 z3 z5 z7) (or x z8))))

2010年4月26日月曜日

HHKBのアニメーション

PFUのHHKB(Happy Hacking Keyboard)のウェブページにはもと和田先生関連ページというのがあって, 私がPostScriptで描いたHHKBの絵が掲載してあった. それをクリックすると, キートップが次々と押されるアニメーションが見られた. そんなプログラムを書いたのは, 1996年頃のことである. (その絵の書き方について, 2004年の夏のプログラミング・シンポジウムで話した. http://www.iijlab.net/~ew/hhkeyboard2.pdf参照.)

最近Processingで動画を描くようになり, その内HHKBアニメも再挑戦しよとう思っていた. 昔のデータがあったので, 実は案外簡単であった.

要するに, 全体の絵と, あるキーが押し込まれたときの, 局所的な絵を用意し, 全体の絵に局所的な絵を重ねるだけである.

全体の絵

Aのキーを押した絵


いま出来ている動画が, http://playground.iijlab.net/~ew/hhkb/hhkb.htmlにあるのでご覧頂きたい. happy hackingと静々と10行打つだけである. 毎秒2文字の早さなので, 70秒で終わる. ただし読み出すのにどういうわけか, いつもかなり時間がかかる. (最新版ではもう少し高速打鍵になっている.)

2010年4月15日木曜日

ケンペの連鎖機構

矢野健太郎:「角の三等分」(ちくま学芸文庫)の注に, 一松先生がケンペの連鎖機構というものがあると書いておられる.

googleで探すと, アニメーションもあるようだが, 私としては, 自分でも描いてみずにはいられない.



Kempeのリンケージは, この図のように, 短い2辺 x と長い2辺 y として, 4本の棒を蝶番でつなげたものらしい. 平面内でくねくねと変形するわけだ. 幾何学的には,角θは0からπまで任意の角がとれる. 実際に物理的なリンケージを作ると, 引っ掛かったりして難しいに違いない.

点Oの座標は(0,0), 点Aは(x,0)である. 点BとCの座標を知りたい. ∠AOB=θとすると, Bは(y cos(θ),y sin(θ))である. すると角φ=tan-1(y sin(θ)/(y cos(θ)-x))と得られる. 従ってBからCを見込む角は2φとなり, Cの座標が得られる. (y cos(θ)+x cos(2φ), y sin(θ)+x sin(2φ)). 座標が得られたら, O-A-C-B-Oと4辺形を描く.

点Cは点Aを中心とした円を描く. 点Bは点Oを中心とした円を描く.その位置を, &thetaの関数と見ると, 点Bは普通の三角関数 cos と sin だが, 点Cの方は, いかにもリンクという曲線になる.



基本形が出来たら, 次にOを中心とし, 座標をθだけ回転し, xの長さを前のyの長さにして, 4辺形を描き, これをもう一度繰り返すと下の図になる. 実はyx の2の立方根倍にして描いてあるので, 最初の x に対し, 最後の y は2倍になっている. 赤い線はそれを示す.



一松先生によれば, このリンケージは, 角の三等分も出来るし, 2の立方根の比もあるので, 2倍の体積の立方体も作れる優れ物らしい.

さて, このリンケージのアニメーションもProcessingで書いた.

http://playground.iijlab.net/~ew/kempe/kempe.html

に置いてある. 上の図のB''に対応する辺りをマウスで動かすと, リンケージ全体が動く.

追記: 教養学部2年の冬学期に矢野健太郎先生(1993年没)の数学の講義があった. その学期のもう1つの数学の講義は彌永昌吉先生(2006年没)の担当であった. 一松信先生はご健在で, プログラミング・シンポジウムに毎年参加されている.

2010年4月9日金曜日

千鳥格子

千鳥格子のコートを着た人を見かけ, 組木細工を思い出した. 数年前に飛騨高山で求めたものだ.



Escherの不思議な絵ではないが, どう出来ているか考えたくなる. 実は予想通りで,


のように, 通常なら幅の1/2だけ削れているほぞが2/3削れている. 最後の1本は1/3だけ削れている(一番下).

図で描くと



上から見たのが 図aである. 交互に上下になって外れそうもないように見える. それぞれの棒にAからFまでの名前をつける. A,B,Cをの下方から見たのが図bで, D,E,Fを右方から見たのが図cである. (第三角法風に描いてある.) 図cに見るようにEはほぞが1/3で, 他の棒は2/3である.

机の上にBを置く. そのほぞに合うように, DとFを載せる. さらにその上にAとCを載せる. 図aに見るように, 確かにA, CはD, Fの上にあり, D, FはBの上にある. そのA, Cの上にEを載せ, それが一番下のBの下をくぐれるであろうか.

この状況で, 一番下のBをつり上げてみる. すると図eで分かるように, F(とD)は, Bより1/3だけ沈む. さらにAとCは, F(とD)より1/3だけ沈む. 従って図dで分かるように, BとCは2/3だけずれ, その間に幅1の穴があくのである. そこにEを滑り込ませればよい.

私の持っている千鳥格子は, このPostscriptの図ほどは正確ではなく, 2/3より多く削れている. だから最後にEを通すのは結構簡単である.

もっと本数が増えても, 理屈は同じだ. 図aで, A, C, D, Fを基数番目, B, Eを偶数番目とし, まず横の偶数番目を置き, 次にその上に縦の奇数番目を載せ, さらにその上に横の奇数番目を載せ, 横の偶数番目を一斉に持ち上げれば, 縦の偶数番目を差し込む穴が出来る. そこへ縦の偶数番目の棒を差し入れればよい.

2010年4月4日日曜日

昼間の長さ

今朝の朝日新聞の「しつもんドラえもん」は夏至と冬至の昼間の長さはどのくらい違うかであった. 14時間35分と9時間45分で約5時間が答. 足して24時間20分になるのは, 日出, 日入の時刻が太陽の上縁が地平線に接する時なのと, 大気差のせいである. 10分引いたくらいが丁度いいかもしれない. ドラえもんの答に北に行くほど差が大きいとあったので, どう変化するか調べてみた. 前のブログ 不等時法にあった式を書き直したのがこれだ.

(define (foo lat) (define delta 23.5)
(let* ((a (sin (d2r delta))) (b (* a (tan (d2r lat)))))
(r2d (asin (/ b (cos (d2r delta)))))))

(define (bar lat)
(/ (* (+ 180 (* 2 (foo lat))) 24) 360))

fooは前のと同じで, θを計算するから, θを昼間の時間に換算するのがbarである.

北緯0度から66.5度までの昼間の時間をプロットしたのが下の図である. 北緯66.5度の北極圏では, 日が沈まなくなるので, 昼間は24時間になる. 赤道上はいつも12時間であるが, こういう傾向とは知らなかった.

またインターネットにある国立天文台の各地のこよみから6月21日の日出と日入の時刻と上の式の結果を書いてみた.

地名 根室 東京 那覇
北緯 43.333 35.658 26.217
日出 3:37 4:25 5:37
日入 19:02 19:00 19:25
昼時間 15:25 14:35 13:48
計算値 15.229 14.423 13.648
時:分 15:39 14:25 13:14

東京は恐ろしいほどよく合っているが, 根室, 那覇は全然違う. 天文台の市川さんに聞いてみようか.

2010年4月2日金曜日

ビットスワップ

2kswapを使い, 64ビットj(j5j4j3j2j1j0)2をjπ=(j0j1j2j3j4j5)2に置換したいというfftのビット変換の解も見事であるが, TAOCPの演習問題7.1.3-53にはもっと驚くことが書いてある.

jπ=(j0j1j2j3j4j5)2をさらに一般的にし, jπ=(j(d-1)ψ...jj)2と書く. ここで{0,1,...,d-1}の置換ψがt個の循環を持ち, d-t個の交換(u1v1)...(ud-tvd-t)について(ただしvi>uiとする), この順の積がψになるなら, この順に&mud,v&¬&muuのマスクで(2v-2u)swapをすればビットスワップが出来るというのだ.

完全シャッフルでは 5ψ=0, 4ψ=5, 3ψ=4, 2ψ=3, 1ψ=2, 0ψ=1だから(0 5 4 3 2 1) -> (5 4 3 2 1 0)には0→1→2→3→4→5→0と1個の循環しかないので, (0 1) (0 2) (0 3) (0 4) (0 5)の交換が必要である. 一方fftでは, 5→0→5, 4→1→4, 3→2→3の3個の循環があり, (u v) = (0 5) (1 4) (2 3)でビットスワップが出来ることになる.

それっというので, 早速前回のスワップをやってみる.
δは25-20=31, 24-21=14, 23-22=4で, マスクは

0000000000000000000000000000000010101010101010101010101010101010

0000000000000000110011001100110000000000000000001100110011001100

0000000011110000000000001111000000000000111100000000000011110000

である.

前回のような図も直ぐ書ける.



前回と殆んど同じ推論が出来る.

31swapでは,
右端の"/"を見るとj0=0, j5=0の時, j0は変らず, j5=0も変らず.
左端の"A"を見るとj0=1, j5=1の時, j0は変らず,j5=0も変らず.
右端の2番目"+"を見るとj0=1, j5=0の時, j0は変る, j5=0も変る.
左端の2番目"B"を見るとj0=0, j5=1の時, j0は変る, j5=0も変る.

これを表にしてみると:



なんのことはない, j0'=j5; j5'=j0であった. つまりビットj0とj5が交換できた. 後は推して知るべし. ビットは左右反転できたのである.

面白いなぁ.

2010年4月1日木曜日

ビットスワップ

3月12日のブログ, ビットスワップの続きである.

2kswapを使い, 64ビットの各ビットj(j5j4j3j2j1j0)2をjπ=(j0j1j2j3j4j5)2で置換したい. そのためのマスクθkを見つけよという問題である(演習問題7.1.3-52).

始めての事ゆえ, 見当もつかず, これもTAOCPの解答を見てみると,

d)0<=k<=5について θk6,k5-k; 0<=k<=2についてθ'kk; θ'3=θ'4=0.

とある. なぜか. 解明したい.

これだけ眺めていても埒があかないので, 64の各ビットがそれぞれのδswapでどこへ移動するかを描いてみた. 幾何の証明の時と同じで出来るだけ正確な図を描くのが肝心である. 幸いとPostScriptがあれば, 一発だ.



この図は
http://www.iijlab.net/~ew/images/bitswapchart.gif

にあるので, これが小さければ, 直接それを見た方がよいかも知れぬ.

δ=1,2,4,8,16,32についての移動パターンは規則的である. 置換πも規則的なので, そうなっても不思議はない.

少し調べてみると, 1swapでは
0****00****1 になり,
0****10****0 になり,
1***** は不変
であることが分かる. *はドントケアを表す.

この伝で他のswapも見たところ

32swap c****t
16swap *c**t*
8swap **ct**
4swap **tc**
2swap *t**c*
1swap t****c
であることが判明した. この読み方で, *はドントケア, tはテストビット, cはチェンジビットで, tが1なら不変, tが0ならcの0と1を変えるのである.

図も規則的だが, ビット変更も規則的である.

さてこれで j5j4j3j2j1j0が j0j1j2j3j4j5になるであろうか.

swapは図のようにδ=1,2,4,8,16,32,4,2,1の順に行う.

1swapではj5=0ならj0は反転し, j5=0ならj0は反転しないから, j'0←j0⊕j5⊕1 と書ける. (変化後のビットにはダッシュ(','')をつけて区別しよう. このように場合を区別するのに, ifを使わず, ⊕などで単一の式で表すのが味噌である.)

同様に, 続く2swap,4swapでj'1←j1⊕j4⊕1; j'2←j2⊕j3⊕1 となる.

次の8swapでも同様にj3が変化するが, これは新しいj'2がtビットなので,
j'3←j3⊕j'2⊕1=j3⊕(j2⊕j3⊕1)⊕1=j2.

同様にして16swapでj'4←j1; 32swapでj'5←j0 となる.

次の4swapではj'2=j2⊕j3⊕1 が新しいj'3つまりj2に従って変化する.
j''2←j'2⊕j'3⊕1=(j2⊕j3⊕1)⊕j2⊕1=j3.

同様にして2swapでj''1←j4; 1swapでj''0←j5 となり, 無事にビットが反転された. QED. めでたしめでたし.