2010年12月26日日曜日

進んで戻って

いよいよ年末が迫る. 元日からの通日は360+αまで来た.

ずいぶん前になるが, 365歩のマーチという曲があり, その中に「3歩進んで2歩さがる」という文句があった. これを聞いて, 昼間に2メートル登り, 夜間に1メートル滑り落ちるカタツムリが, 地上から5メートルの電柱を登りきるのに何日かかるかというクイズを思い出す向きもあろう.

TAOCP V4F3に
"To take one step forward, take two steps forward, then one step backward; to take two steps forward, take one step forward, then another"
という禅問答のような記載があり, 私も直ぐには悟りえず, しばし悩んだことがある.

順列組合せと一緒にしていわれることが多いが, 組合せは, 0,1,...,n-1のn人から, s人を選ぶ方法のことである. これはn - s = t人を選ばない方法でもあるから, TAOCPでは(s,t)組合せという.

つまり, 0号室と1号室があり, 0号室の人数をs, 1号室の人数をtとする. 0番君からn-1番君のそれぞれが, どちらの号室にいるかを, nビットの二進数で表わす.

当然s+tCs=s+tCt組が作れるわけだ.

その組合せを書き出すのに, 常識的には辞書式順, あるいは昇順に並べる. 次は(3,3)組合せの例である.




この図の赤点は, 前後の組合せで異なる位置を示す. 二進数の位置を右から0,1,2,...とすると, 始めの方の最初の赤点は, 000111から001011になった時に, 2番君が1号室から0号室に, 3番君が0号室から1号室に移ったことを示す. 右下の最後の組合せの下に6つの赤点があるのは, 111000から最初の組合せ000111に戻るには, 全員が部屋を移る必要があるということである.

かくも大勢が右往左往しては, ドア回りが混雑するので, 回転ドア方式を考えた人がいる. 0号室と1号室の間に回転ドアがあり, 回転ドアを通って1人ずつだけが入れ替わるのである.



すなわち, この図で見るように, 赤点は常に2個である. なかなかうまく出来ている.

これは6ビットのGray codeを作り, 1のビットが3個のものだけリストして作った.

(define (gray bs)
(define (xor a b) (modulo (+ a b) 2))
(map xor bs (cons 0 (reverse (cdr (reverse bs))))))

(for-each (lambda (n) (for-each display n) (newline))
(filter (lambda (x) (= (apply + x) 3))
(map gray (map (lambda (n) (n2bs n 6)) (a2b 0 64)))))

ところでよく見ると, この赤点の間隔はいろいろだ. 2つ隣り同士のところが10, 間が1個空いているのが4, 2個空いているのが4, 3個空いているのが2で計20である.

この間の幅を0と1にしても出来るといい出した人がいるらしい.

がその1例である. 最後から最初に戻るのは別として, 途中の赤点はすべて隣り同士か1個空きである. このパターンをA33という. しかし, こういうパターンはこれだけでなく, 次のB33もそうなっている.


途中には目をつぶると, A33は111000から始まり, 011100で終る. B33は111000から始まり, 001110で終る. 一般にAstは, 0がs個, 1がt個から始まり, 0が1個, 1がt個, 0がs-1個で終り, Bstは, 0がs個, 1がt個から始まり, 0が2個, 1がt個, 0がs-2個で終る.

AとBの漸化式も知られている.

Ast = 1Bs(t-1), 0A(s-1)tR;
Bst = 1As(t-1), 0A(s-1)t.

これを図で示すと次のようになる. いまはs=t=3である.



左の組合せは, は下の黒文字の説明にあるように, B33である. 上の漸化式に当てはめると, B33 = 1A32, 0A23.

つまり, 左のB33は, A32の前に1を置いたものを並べ, 次にA23の前に0を置いたものを並べると読む. 1や0の後に置く部品の組合せは赤で囲んである. 説明も赤文字である.

また, A33 = 1B32, 0A23R.

つまり, 右のA33は, B32の前に1を置いたものを並べ, 次にA23の, 右肩にRがついているから, 上下逆転したものの前に0を置いたものを並べると読む.

たしかに左の下半分のA23と右の下半分のA23Rは, 上下が逆である.

改めて右のA33を見ると, 最初の111000は最後は011100になり, 111は右へ1歩進んだ. B33では, 111は右に2歩進んだ. 部品を見るとA32は11が右へ1歩進み, A23は111が右へ1歩進み, B32は11が右へ2歩進む. 右下のA23Rを上から下へ眺めると, 今度は01110が11100になるから, 111が1歩戻ったことになる.

これが最初の「一歩進む(Ast)には, 二歩進み(Bs(t-1)), 一歩戻る(A(s-1)tR). 二歩進む(Bst)には, 一歩進み(As(t-1)), もう一歩進む(A(s-1)t)」の謎解きであった.

私がこの解釈に気づいたのは, 2008年の1月頃であった, 2008年の夏にKnuth先生に会った時, 「あの2歩前進, 1歩後退の意味は, 苦労したがついに分かった.」と報告した. Knuth先生は笑っただけであった.

2010年12月25日土曜日

微分解析機

今回は微分解析機といっても, トルクアンプ, 日本語では回転力増幅器の, それも図の話である.

微分解析機の積分機の出力は, 摩擦で生じる微少な力なので, これで他のものを駆動することは出来ない. ところがNewmanにより, トルクアンプが発明され, それによって微分解析機は実用になった.

トルクアンプは, 下の図のような構造である.



右のinputの軸の回転を, 左のoutputの軸から力を増幅して取り出すのである. その間にロープを書けたドラムが2組あり, ロープの端は, 入力と出力軸に取りつけたT字状のレバーの先に固定されている. ドラムは両端のプーリーにかけたベルトで, 矢印の方向に, 逆方向に回転している. 入力軸が回転しない時は, ロープはドラムの面を滑っている.

今, 入力軸が, 右に回転したとする. そうすると, 下のT字状のレバーが持ち上がり, 左のロープがドラムの上で締まり, ドラムとの間に摩擦が生じ, 上のT字状のレバーが押し下げられ, 出力軸が入力軸と同じ方向に回転する. 回転角が同じになれば, ロープは弛み, 再び滑り出す.

入力軸が反対に回る場合は, 右のドラムの摩擦が生じ, 出力軸はやはり入力軸と同回転する.

つまり, 回転するドラムが動力になり, 入力の微少な力を増幅するのである.

こういう仕掛けは船の碇を巻き上げる装置にもあるらしい.

さて, 上のような図で説明したが, 私としては, ロープを直線で描いたのが気に入らないのである. 図はなるべく正確に描きたいというからには, このようにドラムの縁で光の反射のような図はだめである.



そこで考えたのがこの図である. 上と下のT字状のレバーにロープが繋がっている点をAとEとする. 図の円はドラムの断面である. ロープは, AからBに至ってドラムに接し, BCDFBCDとドラムを1回とすこし回り, Dから離れてEでもう一方のT字状のレバーに辿り着く. 先ほどの図は, これを左から見たものだ. そうすると, ロープ上の点は, Aから下がり, ドラムの円の右半分を上がり, また下がるように見えるはずである. また左右の座標位置は, ロープの長さに比例して右へ移動するはずである.

その移動の様子を, 半径の右の, 点の列で示す. 赤はロープのAからBに対応する. 緑はBCDFBCDに対応, 最後の青はDからEである.

そういう解釈で, 最初の図を描き直すと下の絵になる.



ざっと見ると最初の図と左程違わないが, ドラム回りのロープは本当らしくなったのではないか.

2010年12月23日木曜日

クリスマスツリー講義

今月(2010年12月)始め, Knuth先生から, いつもの航空便ではなく, 珍しくメイルで返事が来た. その最後に

By the way, my lecture next Monday will probably be webcast live, so you might be able to watch it via Internet in Japan.

と書いてあった. その日(日本時間では7日朝)は所用あって, インターネットも見ずに過ぎたが, ホームページには,

Monday, 06 December, 5:30pm, in the new NVIDIA auditorium (Huang Engineering Center) A Computer Musing entitled ``Why pi?'' [the sixteenth annual Christmas Tree Lecture]

あ, これがあのクリスマスツリー講義(Christmas Tree Lecture)だったのか.

すぐに連想するのは, 元祖「クリスマス講演」である.

「これから皆さんにロウソクのお話をいたそうと思います.」で始まるMichael Faradayの「ロウソクの科学」は, 1860年のクリスマスであった. だから150年前の話だ. ところでこの原題はThe Chemical History of a Candleという.

Royal Institution Christmas LecturesをWikipediaでみると, 1825年に始まり, 第2次世界大戦中を除き, 今年も続けているというから, 英国のこだわりにはまったく脱帽する.

一方, クリスマスツリー講義は, Knuth先生が, 年末に一般向けに講義するものらしい.

TAOCP V4F4の表4に8次のクリスマスツリーパターンという図があり, その脚注に(TAOCPに脚注があるのは異例なことだが),

"This name was chosen for sentimental reasons, because the pattern has a general shape not unlike that of a festive tree, and because it was the subject of the author's ninth annual "Christman Tree Lecture" at Stanford University in December 2002."

と書いてある.

8次は大きすぎるから遠慮し, 6次のクリスマスツリーは以下のように, 6ビットの二進数64個を, 6C3=20行, 6+1列に配置する.

左からの各列は, 1の個数が0, 1, ..., 6である. 従って中央の高さは6C3である. 同じ行では, 右に行くに従い, 0が1に変っていく. 赤い0は, 次に1になるもの, 青い1は, 前に0だったものだ.

この配置法は, 1951年にde Bruijn, van Ebbenhorst Tengbergen, Kruyswijkが発見した.



101010
101000 101001 101011
101100
100100 100101 101101
100010 100110 101110
100000 100001 100011 100111 101111
110010
110000 110001 110011
110100
010100 010101 110101
010010 010110 110110
010000 010001 010011 010111 110111
111000
011000 011001 111001
001010 011010 111010
001000 001001 001011 011011 111011
001100 011100 111100
000100 000101 001101 011101 111101
000010 000110 001110 011110 111110
000000 000001 000011 000111 001111 011111 111111


これはどうできているかというと, 1次は

0 1

の1行. これではクリスマスツリーにならない. 2次は

10
00 10 11

の2行. いささかツリーらしくなる.

さて, n次のツリーの各行を,
σ1,...,σs
とすると, n+1次のツリーは, 各行を
σ20,...,σs0と
σ10,σ11,...,σs1の
2行にする.

すなわち, 後に0をつけたものと, 後に1をつけたものを作り, 0をつけたものの先頭の1個を1をつけたものの先頭に移動するのである. 先頭を外した後が0個になったら, その行は除く.

Schemeで書くと

(define (next ct)
(if (null? ct) '()
(let ((ss0 (map (lambda (s) (string-append s "0"))
        (car ct)))
(ss1 (map (lambda (s) (string-append s "1"))
        (car ct))))
(set! ss1 (cons (car ss0) ss1))
(set! ss0 (cdr ss0))
(if (null? ss0) (cons ss1 (next (cdr ct)))
(cons ss0 (cons ss1 (next (cdr ct))))))))

である. ctを

(define ct '(("0" "1")))

によって1次のクリスマスツリーとし, ctにnextを順に作用させると,

(("10") ("00" "01" "11"))

(("100" "101") ("010" "110") ("000" "001" "011" "111"))

(("1010") ("1000" "1001" "1011") ("1100")
 ("0100" "0101" "1101") ("0010" "0110" "1110")
("0000" "0001" "0011" "0111" "1111"))


が得られる. クリスマスツリーらしいのは, 次数nが偶数の時である.

6次のクリスマスツリーパターンは, 異なる64の要素があるから, それを八卦の64のパターンで置き換えたのか以下である.



これは, 昨年の, 私からKnuth先生へのクリスマスカードであった.

2010年12月22日水曜日

1のかたまり

TAOCPに  Ank=0n n-kC2k  という式があった(演習問題7.1.4--125).

その解答のHistorical Noteに Ancounts binary x1...xn-1 with each 1 next to another と書いてある. なんだろうと思っているうちに, Knuth先生からのある手紙にヒントがあった.

ちなみにAnの値は

n 1 2 3 4 5 6 7 8 9 ...
An 1 1 2 4 7 12 21 37 65 ... .

とりあえずn=5を考える. するとbinary x1...xn-1

0 0 0 0 *
0 0 0 1
0 0 1 0
0 0 1 1 *
0 1 0 0
0 1 0 1
0 1 1 0 *
0 1 1 1 *
1 0 0 0
1 0 0 1
1 0 1 0
1 0 1 1
1 1 0 0 *
1 1 0 1
1 1 1 0 *
1 1 1 1 *


with each 1 next to another は, 1があれば単独では現れないということで, そういうものは, 上で星印をつけた7=A5個になる. 1は孤立していて失格である. 問題はこのパターンの見つけかたである. 星印のものを取り出すと

0 0 0 0
0 0 1 1
0 1 1 0
0 1 1 1
1 1 0 0
1 1 1 0
1 1 1 1

最初のは1がないから別にして, 残りは1のかたまりが1個ある. 1が孤立するかたまりはないから, 1のかたまりから1個外したパターンを書くと

1 2 3 4
0 0 1
0 1 0
0 1 1
1 0 0
1 1 0
1 1 1

となる. ビットの境界に左から1,2,3,4と番号をつけ, 0と1の境界の位置を書いてみると

3 4
2 3
2 4
1 2
1 3
1 4

となり, これは4個から2個選ぶ数であった. つまり4C2=n-1C2=6.


同様にして, 1のかたまりが2個あるパターンは, n=7の0 1 1 0 1 1 から考え始めて

1 2 3 4 5
0 1 0 1 2 3 4 5
1 0 0 1 1 2 4 5
1 0 1 0 1 2 3 4
1 0 1 1 1 2 3 5
1 1 0 1 1 3 4 5

だから, 右にあるような境界の位置を見れば, 5C4=n-2C4(1のかたまりが2個の場合).

という次第で, 一般に1のかたまりがk個ある数はn-kC2k.
一番上の1のかたまりのないものはnC0=1.



2項係数nCmは, n<mの時に0だから, かたまりがとれるところまで足せばよく, 1のかたまりの置き方の数になるのであった. 結構むずかしい. 私もヒントなしでは, ここまで辿りつかなかっただろう.

2010年12月11日土曜日

2011

1年ほど前 2010が1から9までを使ったどういう式で表わせるかを計算し, このブログに書いたことがある.

今回は2011である. 28個の解があった.


(+ 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)))))
(* 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
が美しい.

2010年11月20日土曜日

双曲線

円錐曲線の中で, もっとも馴染みがないのが双曲線である. したがってda Vinciも双曲線コンパスは考えなかったかも知れない.

反比例 y=1/xを習うと, この形は双曲線と教わる. 確かにx軸とy軸は漸近線だし, 第1象限と同じ図形が第3象限にもある.

xy=1を x'=√2 x-√2 y, y'=√2 7+√2で置き換える. つまり45度回転すると, 確かにx2/2-y2/2=1となり, 双曲線である.

ところで, 双曲線というと, 以前から気になる文言がある. 高橋秀俊先生の「数理と現象」に

「楕円は丸い茶碗やお盆を横から見たときの形として, 万人の日常接しているところであるが, 放物線と双曲線はそれほど親しみがあるとはいえないかもしれない. しかし, 円錐曲線はすべて円錐の切り口, つまり円形をある方向から見たときの形であって, 誰でも目にふれるはずのものである. すなわち大きい円を近くから見て, 円の一部が眼よりも後に来るような場合は, 円は双曲線(の片側)に見え, ちょうど眼の平面(眼のところを通って視線に垂直な平面)に接する場合は放物線に見えるのである.」

と書いてある. 確かに付随する写真では放物線や双曲線のように見える. これはどういうことか, ずーっと疑問に思っていた.

私は自他ともに認めるPostScript描画大好き人間なので, さっそく図を実際に書いてみることにした.




すなわち, 円Oの中央の真上の点Hから, 水平に置いた円Oを, その円周上で接している面Pに投影し, 得られる(赤線の)図形を描いてみたのである. 特殊な環境だが, 計算式を簡単にしたかったからである.

たしかにHの真横の円の部分は右下, 左下の遠方に投影されるから, 双曲線になりそうである. そして出来た図が下である.






視点HをOの真上にとったのは, 怪我の功名であって, 実はHを頂点とする円錐を平面Pで切ったので, 断面が双曲線なのは当たり前である. 高橋先生がいわれるのは, 円内のもっと一般的なところで見ても双曲線ということで, それには別の視点が必要そうだ.

2010年11月17日水曜日

楕円コンパス(つづき)

Leonardo da Vinciは放物線のためのコンパスも考案した. この方はたまたまウェブに図があった(http://www.japandesign.ne.jp/HTM/REPORT/art_review/38/big/05.html)ので, 借用させていただく.



前回のブログの最後に書いたように, 離心率が大きい(=1に近い)と楕円コンパスは苦しくなる. 放物線は離心率が1なので, 発想の転換が必要だ. そこでda Vinciは別の機構を考えたのであろう.

こちも原理は楕円の場合と同じである. 円錐を母線と平行な面で切ると放物線になる性質を利用している. 前回の図で, 切断面が母線と平行になると, 上の球は作図出来るが, 下の球は存在しない. つまり焦点は1つになる.



放物線コンパスの仕掛けを上の図で説明しよう.

ACは鉛直な回転軸である. Cを頂点とし, 底面が円DQEFの円錐を考える. 回転軸を支えるのに, CD, CE, CFの3本がある. 3本の足は, 別に母線に合っている必要はないが, da Vinciの図がなんとなくそうなっているのに従った.

P点に筆記具をつけた伸縮自在の腕 BPは, 母線に沿って回転し, PはDからスタートしてPへ至り, さらにEまで進む. すると放物線DPEが描けるのである.

円錐の角度を一定に保つために, GHがある. 模型で見るように, 放物線を描く面は, 母線と平行に傾けてある. その面の上をPが伸縮しながら放物線を描くのである. 筆記具を平面に押し付けるためか, CPの中間に錘りがぶら下がっているのも面白い.

2010年11月1日月曜日

楕円コンパス

この世には円(お金ではない)も多いが, 楕円はもっと多いに違いない. 丸いものも正面以外の方向から見ればみな楕円に見える. そもそもが, 惑星の軌道から楕円である. また円は楕円の特殊な場合に過ぎぬ.

正弦や余弦の曲線も正確に手書き出来る人はまずいないが, 楕円や放物線も正確には書きにくい曲線である. 最近楕円型の炬燵があると聞いたが, 見ると陸上競技のトラックの形(小判形)であった. 世の常識はそれも楕円だ.

小学生の頃習った楕円の書き方は, 2本のピンと糸を使うものであった. 実際にやってみると, 鉛筆の芯から糸が外れたりで, 結構苦心する. 与えられた長径2a, 短径2bの楕円を描こうとすると, まずピンの間隔を2√(a2-b2), 糸の輪の長さを2(a+√(a2-b2))にしなければならず, 糸の輪の長さをこのように正確にするのは, 極めて難しい.

従って, 製図では, 楕円上の点を適当にとり, 雲形定規でつなげて描いたりする.

最近はPostScriptさまさまで, 縦横のスケールを変えて円を描けば簡単に出来るようになった. しかし線の太さも影響を受けるのが問題である. 長径a, 短径bの楕円は,

a 0 moveto
1 1 360 {dup cos a mul exch sin b mul lineto} for stroke

で描くのがよい. Processingはellipse関数の方が基本で, 円はその特殊なケースとする. この方が正しいやり方である.

しかし, やはりコンパスで円を描くような, 名案はないであろうか.

Archimedesの楕円コンパスというものがある. 我が家にはそのおもちゃがある.



おもちゃとしては, これは何の軌跡かというクイズなのだが. ...

Archimedesの楕円コンパス下の図のようなものだ.



長さaの棒の一端Pに筆記具, 他端と, 筆記具からbの距離に, それぞれ短軸, 長軸に沿って動くスライダN, Mをつけ, その棒を回転すると, 長径2a, 短径2bの楕円が描ける. この動画はここにある.

∠PMQをθ. Pの座標をx, yとする. x=a cos (θ), y=a sin (θ) - (a - b) sin (θ) = b sin (θ)で, Pの軌跡は, 長径2a, 短径2bの楕円になる.

2007年に東京国立博物館で,「レオナルド・ダ・ヴィンチ--天才の実像」という特別展があった. 殆んどの入場者には, 「受胎告知」が関心の的であったろう.

私の場合は, 「レオナルドの書斎」が面白く, 中でもアトランティコ手稿に基づく楕円作図のためのコンパス[複製]に目がとまった.



上の図のようなものである. 通常のコンパス同様, 回転軸ABがあり, Aを固定し, Bをつまんで回す. 通常のコンパスとは異なり, 軸は傾斜している. その傾斜角を維持するために, 支柱CDとCEがある.

軸上の点Qから腕QPが分かれており, Pに筆記具を着ける. AQPの角 θ が一定になるようにして, 軸を回すと, Pが楕円を描くという仕掛けである. ただし, PQは, 軸が傾斜しているから, 通常のコンパスと違い, 伸縮できなければならない.

このコンパスに味噌は, θを一定にして軸を回転すると, QPを母線として円錐が描け, それを平面で切っているのである.

なるほど, 円錐を斜めに切ると楕円が出来ることの応用と思い, さっそくアニメーションのプログラムを書いた.

これは面白かったので, 2007年のlightweight languages spiritで話をした.

円錐を斜めに切ると, 卵形ではなく, 楕円になるというのも, 直ぐには信じ難い.


上のような図を描くといちおう納得する. 円錐を右上から左下へ平面で切った断面である. 円錐と切った平面に上から内接する球を O , 下から内接する球を O'とする. O と円錐の接円, O'と円錐の接円は赤い横線のようになる. また, O , O' と切った平面の接点をそれぞれ F, F' とする.

さて, 任意の母線ABを考える. この母線と Oの接円, 切った平面, O'の接円との交点を, T, P, T'とする. 球面外の点から, 球面へ引いた接線の, 接点までの距離は等しいから, PT=PFであり, PT'=PF'である. 従って, FP+F'P=PT+PT'=TT'=一定となって, Pは楕円となる.

ところで, da Vinciのコンパスの最大の問題は, 伸縮する腕のことだ. 長い時と短い時の長さの比が大きい腕は作るのが困難である. つまり離心率の大きい楕円は描けない. 私の図でも, 比を2:1にとるのがやっとであった. 要するにこのコンパスは話だけで, 実用には遠いのだ.

2010年10月24日日曜日

Maxwellの面積計

またまた面積計の話で恐縮だ.

RutherfordJournal
には, こういう図



も出ている. 前回(10月15日)のブログにあった図と, 似ているがいささか異なる. これもMaxwellのもう1つの面積計なのだ. その理由はご賢察の通りであるが, 老婆心ながら説明しよう.

前回の図では, 水平軸ABを持つ半球が, 車輪Hの回転と共に回転した. また伸縮棒MSにより, 全球と半球を結ぶ直線と水平軸のなす角が変化した.

今回は, 垂直な回転軸に半球がついており, 図形の周囲を追跡する針の, 回転軸の円周方向の移動に応じて, 半球と全球の仕掛けも回転する. また, 針の, 半径方向の移動は, 前回の図のOMの移動同様に, 回転軸と, 半球と全球を結ぶ直線のなす角を変える.

そうと分かると, 前回は車輪の回転がdxで, 半球と全球の角度がyであったように, 円周方向がdθ, 半径方向がrの極座標で, r dθの極座標での積分をすることになる.

前回のMSは伸縮出来る棒であった. 今回も針の先と, 半球の中心と, 全球の中心が1直線になるような, 伸縮棒でもよいわけだが, 今回は曲線部分を持つ部品が回転し, バネにより, それに接触している半球と全球を結ぶ棒の延長が針先になるように工夫したカム機構になっているのである.

Maxwellは理論的な面にしか関心がなかったらしく, ものは作らず, こういうメカを考えて楽しんでいたらしい.

2010年10月22日金曜日

Maxwellの面積計

面積計の歴史をウェブで探すといろいろ見つかる. それを読むと, 円錐(cone)がどうのこうのと書いてあったりする. 面積計の基本知識として, このブログでは円錐の機能を説明しよう.



この図は, Johannes Oppikoferの面積計といわれるものである. Maxwellのと違い, 台枠は固定されていて, 手前にある針を前後に動かすと, その支持腕が伸び縮みし, また腕を左右に動かすと, 支持腕とその右にごちゃごちゃ見える機構が一緒に左右に動く.

ごちゃごちゃの機構の中央に横倒しになった円錐が見える. 円錐は, 機構が左右に動くと, その幅に従った角度だけ水平の軸の周りに回転する. また円錐の上方に, 立っている文字盤があり, その下に小さい円板があって, 円錐に接している. この円板から上は, 支持棒の前後の動きと連動しており, 針が奥に行くほど円板は円錐の頂点に近いところで円錐に接するようになる.

つまり, 針が左右に dx 動くとき, 円板, つまり文字盤の回転する量 dφ は, 針の奥からの距離 y に dx を掛けたものに比例するわけである.



というわけで, 面積計になるのだが, 針を前後に動かす時, 円板が円錐面を横車を押すように動くのが, Maxwellたちの気に入らなかったようである. たしかに, Maxwellの半球と全球の機構はそうなっていない.

また, 円錐の母線が奥の方では低くなるので, 円板や文字盤の高さを変えることも考えなければならない.

Oppikoferの面積計は, 1827年に考案され, 1836年にパリで制作, 販売されたそうだ.

2010年10月18日月曜日

マッチ棒の立方体

すいぶん前のことだが, 数学セミナーにSystem-5というコラムがあった. 東大の高橋秀俊先生, 森口繁一先生, 学習院大学, のちに東大情報科学科の米田信夫先生, 立教大学の島内剛一先生と私が, 持ち回りでなにか書いていた.

今回のブログは, 私が以前そこに書いたものだ.

オランダ Delft大学のvan der Poel先生が大宮の拙宅に来られた時, おみやげといって, 下の写真のようなマッチ棒で出来た箱を下さった. van der Poel先生によると, デンマークのPeter Naur先生(Backus Naur FormのNaurさん)の問題 「接着剤を使わずマッチ棒だけで立方体を作れ」を考えた結果作ったものだそうだ.

van der Poel先生に頂いたものは, 遠の昔にどこかへいってしまった. この写真のは, 私がこの週末に作ったものである.



マッチ棒を90本も使っているから, ちらちらするが, 三面図を書いてみると, このようになる. 第3角法になっていて, 右下が正面図, 左が側面図, 上が平面図である. マッチの軸に色が着いていて, 薬の部分が白くなっている.




原理的には, 下の図のように, 3対のカードを組合わせ, 対同士を輪ゴムで止めると立体が出来るのと同様である.



実際に作ってみると, この写真のようになる.



ただ, 輪ゴムは, 接着剤ではないが, 違反であって, 輪ゴムの代りにマッチ棒を使えばいいのである. (上の図の左下.) この図ではマッチ棒は抜けてしまいそうだ. しかし, 別の面の力で絞められると, 抜けなくなるのである.

というわけで, 下の図のように, 1,2,3,...,21の順で, マッチ棒を組み立てると, この立方体が完成する. この図では, 薬の方に色がついている. 20と21では, 上に薬が来る時に色がついている.





この図と三面図の正面図の対応は, こうだ. 1は正面図の底の緑の横向きのマッチ. 2は下から2段目の青の奥向きのマッチ. 3から17は, 赤の横向き, 奥向きのマッチ, 18は上から2段目の青の奥向きのマッチ. 19は蓋の緑の横向きのマッチ. 20と21は平面図の奥向きのマッチである.

しかし, 実際やってみると, 賽の河原よろしく, 出来上がる前に崩れてしまう. まぁ, 完成するには, よく考え, ちょっとしたノウハウを利用する必要がある. それは ひ み つ.

2010年10月15日金曜日

Maxwellの面積計

Maxwellといえば電磁気学で習う4つの式や悪魔(demon)を思い出す. 昔の学者はなんでも屋で, そのMaxwellは面積計も, しかも2種類も, 発明したことを最近知った.

Internetで探すと, The Scientific Letters and Papers of James Clerk Maxwellにあるこういう図


や, RutherfordJournalにあるこういう図



が見つかる.

今回は, この面積計の説明である.

下の図を見て欲しい.


これは面積計を上から見たところである. 水平に置かれた回転軸ABがあり, 右方に車輪Hがある. 左方には, 半球KCBがあり, Bが極に対応する. 手前の方, VXに沿って動くMと, 半球の中心Sを通る鉛直軸を結ぶ, 伸縮する棒MSがある. MSから直角に出た枠SDがあり, EFを回転軸とする全球ECB'Fが半球にCで接している. MがOの位置にあると, SDはASと重なり, 全球の赤道DB'が半球の極Bに接し, MがOから左へ行くにつれ, 全球はSを中心に上にあがるが, 球面同士は滑らないものとする. 従って, ∠OSMをθとすると, ∠BSC, ∠B'DCもθである.

面積を計測する時は, 面積計を平行に手前に引いたり, 遠方に押したりしながら, Mが図形の周囲を辿るようにする. RutherfordJournalの図の破線や, 台車の右下の車輪から推察出来る通りである. 面積計の押し引きに応じ, 車輪Hが紙との摩擦で回転し, 同軸の半球も回転し, それに接する全球も回転する.

OMの距離をy, OSの距離をh, 車輪Hの半径をR, 半球, 全球の半径をr とし, 面積計をdxだけ奥に押した時の全球の回転角dψを計算してみる. xの正方向は, Hと半球の向こうへの回転, 全球の手前への回転とする. Hの回転角(radian)は dφ=dx/R. するとC点で半球の緯線が鉛直方向に動く距離は, rsinθdφ. 全球のCでの緯線の半径はr cosθ. 全球の緯線も鉛直方向に同じ長さ動くから, 従って回転角dψ=r sinθdφ/r cosθ=dφtan θ. ところで, tanθ=y/h. ∴dψ=dφy/h=(dx y)/hR.

この面積計で円の面積を計ってみる. R=1, h=2, 円の半径=1とする.



最初, Mを計測対象の赤線の円の真下に置く. この時の半球と全球の中央の経線を赤で示す. 以後の図で,この経線の移動が分かる.



円周の1/4のところまでMが移動した. 面積計も奥へ移動し, 車輪Hが回転し, 半球の目印の経線も1radianだけ回った. ここまででスキャンした面積をオレンジ色で示す. それに相当した角度だけ, 全球も回転した.



Mは円の真上まで来た. 面積計はさらに奥まで移動した. 半球も全球も回転を続ける. 半球の回転角は2 radian. 半球の下側に入ったので, 破線で示す.



円周の右を辿って下がり, 3/4来たところ. 車輪も半球も先ほどとは逆に回転し, 面積を記録している全球も逆回転した.



Mは円周を一周した. 半球の経線は元に戻り, 全球の経線は, 面積πの半分を円周に一致して示す.

Maxwellはこの面積計を作成したわけではない. 実際に作るのは難しそうである.

2010年9月29日水曜日

学者猿コンサル

前回のブログ(コンサル)は, 検討も設計も多少杜撰であった. 改めて考えた描き方はこうだ.

これは2変数関数の値を表示するいわば計算図表みたいなものである. 計算図表には, 変数の変域があるわけで, 今回はそれをxminと, xmaxとする. また2つの変数をx0と, x1とする. つまりxmin≤x0, x1≤xmaxである.



xminとxmaxの基線の上に, x0とx1の点A, Bをとり, それを脚とする長さlの二等辺三角形CABを描く. ACをCで左に135度曲げてDの方向へ長さm(=l/√2)だけ伸ばしDとする. BCをCで右に135度曲げてFの方向へ長さmだけ伸ばしFとする. CD, CFを2辺とする菱形CDEFを作ると, Eが関数の値の場所になり, その座標は(x0+x1)/2, (x1-x0)/2となる.



x座標は誰にでもすぐ分かるが, y座標がかくもきれいな形なのは不思議である. これからは幾何学だが, まず二等辺三角形CABの底角はどちらもαなので, 頂角∠ACB180-2αである. ACの長さはl, ∠ACDは45度, CDはl/(√2)なので, DACは直角二等辺三角形である. DEの線はDから下向きに45-αなので, 菱形で上半分の角度も同じゆえ, ∠CDEは2α-90. 従って, ∠ADEは180-2αとなり, 先程の頂角とおなじである. また, DA=DE=mだから, DAEはCABと1:√2の相似形であった. 従って, AEはABの1/√2であり, 右側も同じなので, EABは直角二等辺三角形である. ゆえにEのy座標, EHはABの半分, つまり(x1-x0)/2である.

y座標は, x0≤x1なら正だが, 反対だと負になり, 基線より下になることに注意しよう.


こうして作った減算表(x1-x0)が次のものだ. 8-3=5と3-8=-5を示している. もうアニメーションを作る必要はないと思う.


2010年9月26日日曜日

学者猿コンサル

2010年9月 Brisbaneで開催されたIFIP WCC2010でHistory of Computingのシンポジウムがあった. そこでオーストラリアのMonash大学のJudy Sheardさんが同大学のコンピュータ博物館を紹介した. そのスライドの中に, おや?と思う絵があった. 一瞬の内に次に進んだが, 同行の山田さんがその絵を写真に撮っておられたので, 頂いた.

Consul, the Educated Monkeyという計算関連の教育玩具である. ウェブページで調べると, 1916年頃にアメリカで売り出されたものらしい.


(出典http://www.rechenwerkzeug.de/consul.htm)

使い方はこうだ. 猿の両足を, 図のように下の目盛の3と9に合わせると, 猿の手が示す枠の中に積の27が現れるのである. 目盛の右端の12の右には四角が見えるが, 片足をそこに合わせると, 他方の足の自乗が現れる.

これは乗算表であるが, この表は差し換えられ, 加算にも使える.

こういうものがあると知ったら, プログラムしたくなるのは当然である.

http://playground.iijlab.net/~ew/consul/consul.htmlを見て欲しい. これは十六進の乗算表である. 下の目盛は0からfまであり, 目盛に乗っている足の辺りをマウスで横に動かすと, リンク機構がつられて動き, 赤印の交点の直ぐ上に積が現れるようになっている.


多少の設計ミスで, 0とfの積などでは, 赤印が頭の上に行ってしまうのは, ご愛敬である.

なお, Consulのシミュレータもウェブにあった.

http://www.edumedia-sciences.com/de/a572-consul-the-educated-monkey

2010年9月16日木曜日

入れ子のかっこ

前回のブログの点記法の続きである. とりあえず練習をしよう. 私の手元のPaul Rosenbloom, The Elements of Mathematical Logic (Dover 1950)に点記法の論理式が沢山出てくる. それでテストしてみた. まずSchemeでBoole値が0と1のnot, or, impを定義する.

(define (not p) (- 1 p))
(define (or p q) (quotient (+ p q 2) 3))
(define (imp p q) (or (not p) q))

テストをするには, mapが便利.

(map not '(0 1)) => (1 0)
(map or '(0 0 1 1) (0 1 0 1)) => (0 1 1 1)
(map imp '(0 0 1 1) '(0 1 0 1)) => (1 1 0 1)


上の本では, 沢山並ぶ -> には特にかっこを(点で)示さない. -> はleft associativeであって, (p -> q -> r) は ((p -> q) -> r) なのだ.

ではやってみよう.

I4 p -> .q -> r. -> .p -> q -> .p -> r
かっこに変えると
((p -> (q -> r)) -> ((p -> q) -> (p -> r)))
Schemeの定義は
(define (i4 p q r)
(imp (imp p (imp q r)) (imp (imp p q) (imp p r))))
実行
(map i4 '(0 0 0 0 1 1 1 1) '(0 0 1 1 0 0 1 1)
'(0 1 0 1 0 1 0 1))
=> (1 1 1 1 1 1 1 1)

もう1つ.

T2 q -> r -> .p -> q -> .p -> r
かっこに変えると
((q -> r) -> ((p -> q) -> (p -> r)))
Schemeの定義は
(define (t2 p q r)
(imp (imp q r) (imp (imp p q) (imp p r))))
実行
(map t2 '(0 0 0 0 1 1 1 1) '(0 0 1 1 0 0 1 1)
'(0 1 0 1 0 1 0 1))
=> (1 1 1 1 1 1 1 1)

とうまくいく.

かっこ記法と点記法の変換だが, まずかっこから点へは, すべての2項演算子をかっこでくくることにし, つまり

<primary>==<letter>|(<expression>)
<expression>==<primary>|<expression>*<expression>

だけとする. expressionの例は a, a * b, (a * b) * c, ...など.

expression * expression を点記法にするには, expression' lp * rp expression とする. expression'はexpressionを点記法に変えたものである. lpとrpは, 左右のexpressionをひとまとまりにする点である. lpは左のexpression'で使った最大の右点の数より多く, rpは右のexpressionで使った最大の左点の数より多くなければならない. 従って, 下請けの変換は, expression'を返すと同時に, 自分の使った右点, 左点を返すことにする.

letterの場合は, letterの他に, 右点,左点として, 0,0を返す.

* の場合は, 左の式を変換て置き, その右点+1の左点を置き, 演算子を置き, 右の式を変換し, その左点+1の右点を置き, 右の式を置く. また自分で使った左点, 右点も返す.

(define (dotconv exp)
(display (list exp))
(if (symbol? exp) (list '(0 0) exp)
(let* ((l (dotconv (car exp)))
(r (dotconv (caddr exp)))
(op (cadr exp))
(rlp (+ (cadar l) 1))
(le (cadr l))
(lrp (+ (cadar r) 1))
(re (cadr r)))
(list (list rlp lrp)(list le rlp op lrp re)))))


(dotconv 'a) => ((0 0) a)
(dotconv '(a * b)) => ((1 1) (a 1 * 1 b))
(dotconv '((a * b) * c)) => ((2 1) ((a 1 * 1 b) 2 * 1 c))
(dotconv '((a * b) * (c * d))) =>
((2 2) ((a 1 * 1 b) 2 * 2 (c 1 * 1 d)))

この後は ((a 1 * 1 b) 2 * 2 (c 1 * 1 d)) を a * b . * . c * d にしたい. flattenし, 整数はそれ引く1の点を出力する.

(flatten '((a 1 * 1 b) 2 * 2 (c 1 * 1 d))) => (a 1 * 1 b 2 * 2 c 1 * 1 d)

(define (convstream l)
(apply string-append (apply append
(map (lambda (x)
(list (if (number? x) (make-string (- x 1) #\.)
(symbol->string x)) " "))
(flatten l)))))

(convstream '((a 1 * 1 b) 2 * 2 (c 1 * 1 d)))

(define (par->dot exp) (convstream (cdr (dotconv exp))))

(par->dot 'a) => "a "
(par->dot '(a * b)) => "a * b "
(par->dot '((a * b) * c)) => "a * b . * c "
(par->dot '((a * b) * (c * d))) => "a * b . * . c * d "


一方, 私の考えた逆変換はこうだ.


1.(((a * b) * (c * d)) * ((e * f) * (g * h)))
2."a * b . * . c * d .. * .. e * f . * . g * h "

1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2
0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8
3.(a 0 * 0 b 1 * 1 c 0 * 0 d 2 * 2 e 0 * 0 f 1 * 1 g 0 * 0 h)
4.29
5.((-1 5) (-1 13) (15 21) (23 29) (15 29) (7 13))
6.(< < < a * b > * < c * d > > * < < e * f > * < g * h > > >)
7.(((a * b) * (c * d)) * ((e * f) * (g * h)))


上の例で, 1. は元のかっこ記法の式. 2. はそれを点記法にしたもの. 逆変換はここから始まる. まず各演算子に左点と右点があるものとし, 3のように変換する. その上の2行は, 各要素の位置を示す. 0から28まであるから, lengthをとると, 4. のように, 29.

位置 1,5,9,..のように, 4を法として1の位置は左点. 3の位置は右点である. さらに, 演算子は2(mod 4), 変数は0(mod 4)である.

次に各右点>0にはこの右方にある相棒の左点を探し, また左点>0にはこの左方にある右点を探し, 対にする. 上の例では 5. が対のリストである. この読み方は, 5の位置の左点を越える右点は, 左方にはないので, -1とし, (-1 5)とする. 13の位置の左点も同じ. 21の位置の左点は, 15の位置の右点の方が大きいので, スコープはここまでとなり, (15 21)が出来る. 3. と5. の情報から, 6. を作るのだが, (-1 5)のような対があれば, -1に左かっこ, 5に右かっこを置く. この挿入で, 位置がずれると困るから, 挿入は番号の多い右端から行う. 右かっこと左かっこの区別は, 4を法とした剰余の1か3で決る. まだこの段階は記号列であるが, 通常のかっこを記号として挿入すると, 分かり難いから, 角かっこを使っている. 出来たのは6. である.後は, 入れ子の式の読込みルーチンを書けばよい. それにより, 7. が得られる.

(define (iconv s) ;2. -> 3.
(define (char->symbol c)
(string->symbol (list->string (list c))))
(define (reads n s)
(cond ((null? s) s)
((char=? (car s) #\Space) (reads n (cdr s)))
((char=? (car s) #\.) (reads (+ n 1) (cdr s)))
(else (cons n
(cons (char->symbol (car s)) (reads 0 (cdr s)))))))
(let ((ss (string->list s)))
(cons (char->symbol (car ss)) (reads 0 (cdr ss)))))

(define (makepair s) ;3. -> 5.
(let ((len (length s)) (ps '()))
(do ((i 3 (+ i 4))) ((>= i len))
(let ((r (list-ref s i)) (k len))
(do ((j (- len 4) (- j 4))) ((< j i))
(let ((l (list-ref s j)))
(if (< r l) (set! k j))))
(if (> (- k i) 2) (set! ps (cons (list i k) ps)))))
(do ((i (- len 4) (- i 4))) ((< i 0))
(let ((l (list-ref s i)) (k -1))
(do ((j 3 (+ j 4))) ((> j i))
(let ((r (list-ref s j)))
(if (> r l) (set! k j))))
(if (> (- i k) 2) (set! ps (cons (list k i) ps)))))
ps))

(define (makestr ps s len) ;5. 3. (length 3) -> 6.
(define (insert l n a)
(if (<= n 0) (cons a l)
(cons (car l) (insert (cdr l) (- n 1) a))))
(set! ps (sort (apply append ps) >))
(for-each (lambda (x) (set! s (insert s x
(if (= (modulo x 4) 3) '< '>)))) ps)
(append '(<) (filter symbol? s) '(>)))

(define (str->sexp str) (define i 0) ;6. -> 7.
(define (getch) (let ((ch (list-ref str i)))
(set! i (+ i 1)) ch))
(define (read) (let ((ch (getch)))
(cond ((eq? ch '<) (readtail))
((eq? ch '>) '())
(else ch))))
(define (readtail) (let ((x (read)))
(if (null? x) x (cons x (readtail)))))
(read))

2010年9月5日日曜日

入れ子のかっこ

Lispの好きな人は, かっこだらけのS式を何とも思わない. 反対にhtmlのようなタグでの入れ子には耐えられない.

しかし, S式を声を出して読む時, 「かっこ開く」「かっこ閉じる」や「左かっこ」「右かっこ」というのは煩わしい. 私は昔から,左かっこを「かっこ」, 右かっこを「こっか」と読むことにしている. (Algol 68で, 添字に使う [ と ] を, sub symbol, bus symbolと呼んでいたのをまねした.)

ところで, 今回の話題は, 自然言語の構文解析である. 構文解析木を図で表わすのは容易だが, かっこ付きでも表わせる.

Time flies like an arrow は 光陰矢の如し と構文解析するのが普通だが, 時間の蝿は矢が好きだ という笑い話もある. 前者は Time (flies (like an arrow)), 後者は (Time flies) like an arrow と表現すればよい.

SICPに簡単な自然言語解析の話がある. その問題4.45は, The professor lectures to the student in the class with the cat を5通りに構文解析せよ, だ. やってみると

1 ((the professor)
(((lectures
(to (the student)))
(in (the class)))
(with (the cat))))
2 ((the professor)
((lectures
(to (the student)))
(in ((the class)
(with (the cat))))))
3 ((the professor)
((lectures
(to ((the student)
(in (the class)))))
(with (the cat))))
4 ((the professor)
(lectures
(to (((the student)
(in (the class)))
(with (the cat))))))
5 ((the professor)
(lectures
(to ((the student)
(in ((the class)
(with (the cat))))))))

それぞれの和訳は

1 教授は(猫を連れて(教室で(学生に講義する)))
2 教授は((猫を連れた教室)で(学生に講義する))
3 教授は(猫を連れて((教室の学生)に講義する))
4 教授は((猫を連れた(教室の学生))に講義する)
5 教授は(((猫を連れた教室)の学生)に講義する)

猫を連れているのは, 1と3は教授. 2と5は教室. 4は学生である.

ウェブを眺めていて, H.B.Curryの書いたページを見つけた. ただで見られるのは1ページ目だけだが, 今はそれだけ見られれば充分である.

論理の式には, かっこの代りに点が打ってあるものがある. 件のページには, その提案が書いてある.

Let us suppose that a group of dots on the right of an operation or prefix denotes the beginning of a bracket which extends to the right until it meets a group with an equal or larger number of dots on the left of an operation; and that the scope of a group of dots on the left of an operation shall extend to the left unitl it reaches a larger group of dots on right of some operation.


つまり,
p -> q . -> . q -> r .. -> . p -> r

((p -> q) -> (q -> r)) -> (p -> r)
のことである.

この方式を何と言うか知らないが, dot notation つまり 点記法ということにする.

上の和訳の入れ子のかっこ構造から, 点記法を使い, かっこを外したい. それには, 演算子がどれか決めたい. 演算子は左かっこの右, 右かっこの左にあるのだから,

教授猫を連れ教室学生講義する

の赤字の助詞類を演算子と考える. (猫を連れた, 教室の は形容詞句の時で, 副詞句の時は, 猫を連れて, 教室で になる.)

すると, それぞれの例文は

1 教授は . 猫を連れて . 教室で . 学生に講義する
2 教授は .. 猫を連れた教室 . で . 学生に講義する
3 教授は . 猫を連れて .. 教室の学生 . に講義する
4 教授は ... 猫を連れた . 教室の学生 .. に講義する
5 教授は .. 猫を連れた教室 . の学生 . に講義する

となる.

この点を , 2点を ,, に変えると
1 教授は, 猫を連れて, 教室で, 学生に講義する
2 教授は,, 猫を連れた教室, で, 学生に講義する
3 教授は, 猫を連れて,, 教室の学生, に講義する
4 教授は,,, 猫を連れた, 教室の学生,, に講義する
5 教授は,, 猫を連れた教室, の学生, に講義する
のようになる.

点記法に慣れると, これで読めなくもない. では声を出して読むにはどうするか.

Shakespeareの時代, 劇の台本には句読点が何種類かあり, それぞれで間の取り方が違うという話を読んだことがある. そのように, カンマ2つはカンマ1つより長く間をとれば, この例は構文木が分かるのではないかと考えるが, どうだろうか.

2010年8月29日日曜日

Penroseタイル

GoogleでPenrose tileの記事を眺めていたら, Penrose tileは3 colorableという記述が! ジャジャーン. 平面の地図は4色で塗り分けられるが, 3色で塗り分けられるというわけだ. さっそく前回のブログの最後の図をプリントし, 4色のマーカーを手に, 塗り分けてみると, 確かに3色でうまく行く.

余談だが, 慶應の皆さんは, 慶應の校旗はblue red and blueの3色といい, 私の2010年5月13日のブログのように, オランダ国旗も3色だ. ジャジャーン.言語学者の用語では, 「サイタ サイタ サクラ ガ サイタ」は, 延べ語数5,異なり語数3という. その伝では, 慶應の校旗は述べ3色旗. オランダ国旗は異なり3色旗である. Penrose tileは異なり3色で塗り分け可能である.

数学者なら, 3色塗り分け可能を証明するわけだが, 私の場合はプログラムで,3色塗り分けしたい. さっそくトライ.

前回のブログの最後にプログラムがあるが, そこで分かるように, それぞれの線分はtranslateで原点を移動し, rotateで向きを変えたユーザ座標て描いている. これでは誰が隣りかは判明しない. どうすれば絶対座標(装置座標)で描けるか, PostScriptのマニュアルをみると, transformという関数がある. PostScriptでは,装置座標とユーザ座標の関係は, current matrixが記憶していて, transformはユーザ座標を装置座標にしてくれるらしい.

ところが,例えば kite0のn=0の部分を0 0 transform == phi2 36 xy transform == phi -72 xy transform ==のようにし,
[1 0 0 1 0 0] setmatrix
300 400 translate
/n 0 def init
0 0 0 kite0

で起動すると, kite0の3点は,

A 300.0 400.0
B 469.442719 523.107361
C 340.0 276.892639
と出力された.

今, 単位長さが80なので, すこし計算してみる.
(* 80 phi phi (cos (d2r 36))) => 169.4427190999916
(* 80 phi phi (sin (d2r 36))) => 123.10734148701015
(* 80 phi (cos (d2r -72)) => 40.00000000000001
(* 80 phi (sin (d2r -72))) => -123.10734148701013

理由は分からないが, AとBは正しい; CはAから-72度の方向に移動した距離になっている. それが分かればあとはプログラムでCの位置は計算出来る.

こうしてkite0, kite1, dart0, dart1の3点の座標が得られた. 次はkite0とkite1, dart0とdart1を組合わせる. 始点と終点は合っているから, それらを対にする. 相棒のないのは外す. こうして95の四辺形が揃った. 各四辺形には0から94まで番号をつける.

次は隣りどうしを探す作業で, alistを作り, 辺を登録したり探索したりして,隣組のリストが出来た. 両方から指すこともないので, 自分より番号の少ない隣りの四辺形のリストになっている.

次にそれぞれに色0,1,2を対応させる. 目で見てやるのと違うから,バックトラックすることになる. SICPの8クィーン(ex2.42)のプログラムをちょっと改造して書いてみたが, これは全解探索だし, 深さも8段から95段になるので,たちまち再帰のオーバーフローを起した. うーん. 駄目か.

仕方ないので, 通常のバックトラックのプログラムを書く.解が1つ見つかったところで, call-ccを使って脱出する.

(define next' (() (0) () (0 2) (2 1) (3) (5) () (7) ()
(7 9 0) (10 0) (9 8) (7) (7 13 5 3) (6 13) (15 6) () (17) ()
(17 19) (19 18) (20) (22) () (24 2) (25 2 5) (4) ()
(24 28 17) (29 17) (28 27 25) (24 26) (24 32 22 20)
(23 6 32) (16 34 23) () (36) () (36 38) (38 37) (39) (41) ()
(43 19) (44 19 22) (21) () (43 47 36) (48 36) (47 46 44)
(43 45) (43 51 41 39) (42 23 51) (35 53 42) () (55) ()
(55 57) (57 56) (58) (60) () (62 38) (63 38 41) (40) ()
(62 66 55) (67 55) (66 65 63) (62 64) (62 70 60 58)
(61 42 70) (54 72 61) () (74) (8) (76 8 13) (74 76) (76 75)
(79 12) (78 77) (81 15) () (83 57) (84 57 60) (59) ()
(83 87 74) (88 74) (87 86 84) (83 85) (83 91 81 78)
(82 61 91) (73 93 82 16)))
(define colorlist '())
(define (search)
(call-with-current-continuation
(lambda (throw)
(define (safe? n c)
(let ((es (list-ref next n)) (revc (reverse colorlist)))
(not (member c (map (lambda (x) (list-ref revc x)) es)))))
(define (test n)
(if (= n 95) (throw (reverse colorlist))
(do ((m 0 (+ m 1))) ((= m 3))
(let ((c (modulo (+ m n) 3)))
(if (safe? n c)
(begin (set! colorlist (cons c colorlist))
(test (+ n 1))
(set! colorlist (cdr colorlist))))))))
(test 0))))
(search)

ジャジャーン. これでやっと色のリストが得られた.

(0 1 2 1 0 2 0 1 2 0 2 1 1 2 0 1 2 2 0 1 0 2 1 2 1 0 1 1 1 0
1 2 0 2 1 0 0 1 2 1 0 2 0 1 2 0 1 1 2 1 0 2 0 1 2 1 2 0 2 1
0 1 2 0 1 2 2 0 2 1 0 1 2 0 2 0 1 0 0 2 0 1 2 0 1 2 2 2 1 0
0 1 2 0 1)

結果はご覧の通り.

2010年8月25日水曜日

Penroseタイル

今回はPenrose tilingの絵の描き方のことである. Penrose tileには

のようにkite(凧)とdart(矢)の2種類のタイルを使うものと,

のようにthick(でぶ)とthis(やせ)の2種類のタイルを使うものとがあるらしい.

これらのタイルは, 菱形以外に組合わせると, 平行移動方向の周期的パターンで平面を埋めることが出来ないことで知られている. どちらのタイルにも, 赤と青の円弧が描いてあるのは, 菱形が出来るのを禁止するからで, タイルを置く時は, 同じ色の円弧が接続する必要がある.

従って, 沢山のタイルによるパターンを描くには, タイルを同型のタイルに分割し(これをdeflationという), 再帰的に埋めていく方法を用いる.

分割はタイルの半分について行うので, 分割を示す破線が描いてあり, 緑の線は半分のタイルの起点を示す.

半分のタイルは, 上には0, 下には1をつけて区別する.

今回は, 私がPostScriptで書いたプログラムの話をしたい.

上の2組のタイルの絵は, googleで探したもので, 凧の方はWolfram MathWorldのだ. kiteとdartが長さで与えられ, thickとthinが角度で与えられているのも, 出所が違うからだ.

kiteの辺ABがφ+1, BCがφ-1+1と示されているが, WolframMathWorldにそう書いてあったからだが, kiteの青の円弧の半径がφで, 赤の円弧の半径が1, dartの青の円弧の半径が1で, 赤の円弧の半径がφ-1であることも示す. φはもちろん黄金比 (√5+1)/2である. φの計算に慣れている人には分かるように, φ+1は実はφ2, φ-1+1は実はφである. 従ってAB:BCはφ:1である. kiteとdartは, dartを左右逆にするとkiteに密着出来て菱形になる. 菱形の長い対角線をφ:1に分けた点がCになる. つまりkiteのAC:dartのAC=φ:1. kiteのCABもdartのBCAも2等辺三角形である.

kiteの∠CABを計算しよう. AからBCの中点Eに垂線を引くと, 直角三角形AEBについて, sin ∠EAB=1/2φ=sin 18°, つまり∠CAB=36°である. ゆえに∠ABC = ∠ACB=72°である. dartについては, ∠CAB=∠ABC=36° である.

(注: sin 3α=3 sin α-4 sin3α, cos 2&alpha=1-2 sin2α. &alpha=18°の時, sin 54°=cos 36°からsin 18°は解ける. sin 18°は(√5+1)/4.)

ここでthinの上半分を見ると, ∠CAB=∠ABC=72°だから, thin0はkite0と同じ形である. また大きさは違うが, thick0はdart0と相似形である.

kite0は下のようである. この図の上は, 1個のkite0の描き方で, 左端の原点(緑の線が出発する位置)から, 36°方向へφ2だけ線を引き, 次にそこから-72°方へφだけ引くことを示す.

図の下は, 1個のkite0の分割法で, kite0'が出来るところである. (分割後のタイルにはプライム(')をつけて区別する.)kite0'の左の三角は, dartの片割れだが, 緑の線からdart1であることが分かる. また右の方はkite0とkite1で出来ている. つまり左端に原点を持ち, 36°に向けたdart1を置き, 36°方向へφ2行った所に原点を持ち, 252°に向けたkite0とkite1を置けばよいことを示す. kite0は長い辺がφ, 短い辺が1であった. 分割の結果は, 短い辺がkiteの長い辺で出来ているから, タイルの大きさは1/φになった. 従って, 分割の図のスケールはφ倍で示してある. 上の図のφ2が下の図ではφ3になるのは, そのためである.

そろそろプログラムの出番だ.

描画の出発点は緑の線のついている頂点で, 頂点の座標x, yと緑の基本の線の方向aを貰い, タイルを描く. kite0は(0,0)から出発, 36度の方向にφの長さの線を引く. 次に-72度の方向に1の長さの線を引く. これをPostScriptで

0 0 moveto φ 36 xy lineto one -72 xy rlineto

のように書く. xyは線の長さ l と方向 a を貰い, その先のx, yの位置を返す関数である.

/xy {2 dict begin /a exch def /l exch def l a cos mul l a sin
mul end} def

kite0'の描き方は, まずスケールを1/φにする. (0,0)の位置に36度の方向でdart1を描く. 次は36度の方向にφ行った所を原点とし, -108度の方向でkite0とkite1を描く. 注意点は, φ行ったところといっても, スケールが1/φになっているから, φ2行く必要があることだ. このことさえ忘れなければ, すべての絵は描ける.

0 0 36 dirt1 phi3 36 xy 252 kite0 phi3 36 xy 252 kite1

phi3はφ3に相当する長さのつもりである.

kite1, dart0, dart1の分割は以下の通り.

dart0

dart1

thickとthinも同じだ.
thick0

thick1

thin0

thin1


kiteとdartのプログラムは以下のとおり.

/phi 5 sqrt 1 add 2 div def /one 80 def
/init {/phs n 5 add array def /f one def
n -1 0 {/i exch def /f f phi div def phs i f put} for
phs n 1 add one put /one one phi mul def phs n 2 add one put
/one one phi mul def phs n 3 add one put
/one one phi mul def phs n 4 add one put} def
/ph {phs exch n add get} def
/xy {2 copy cos mul 3 1 roll sin mul} def
/red {1 0 0 setrgbcolor arc stroke} def
/blue {0 0 1 setrgbcolor arc stroke} def
/green {0 1 0 setrgbcolor 0 0 moveto 1 ph 0 lineto stroke}def

/kite0 {3 dict begin /a exch def /y exch def /x exch def
gsave x y translate a rotate
n 0 eq {0 setgray 0 0 moveto 3 ph 36 xy rlineto
2 ph -72 xy rlineto stroke
0 0 2 ph 0 36 blue 3 ph 0 1 ph 108 180 red green}
{/n n 1 sub def 0 0 36 dart1 4 ph 36 xy 252 kite0
4 ph 36 xy 252 kite1 /n n 1 add def} ifelse
grestore end} def

/kite1 {3 dict begin /a exch def /y exch def /x exch def
gsave x y translate a rotate
n 0 eq {0 setgray 0 0 moveto 3 ph -36 xy rlineto
2 ph 72 xy rlineto stroke
0 0 2 ph -36 0 blue 3 ph 0 1 ph 180 252 red green}
{/n n 1 sub def 0 0 -36 dart0 4 ph -36 xy 108 kite0
4 ph -36 xy 108 kite1 /n n 1 add def} ifelse
grestore end} def

/dart0 {3 dict begin /a exch def /y exch def /x exch def
gsave x y translate a rotate
n 0 eq {0 setgray 0 0 moveto 3 ph 36 xy rlineto
2 ph 252 xy rlineto stroke
0 0 1 ph 0 36 blue 2 ph 0 0 ph 72 180 red green}
{/n n 1 sub def 0 0 0 kite0 4 ph 36 xy 216 dart0
/n n 1 add def} ifelse
grestore end} def

/dart1 {3 dict begin /a exch def /y exch def /x exch def
gsave x y translate a rotate
n 0 eq {0 setgray 0 0 moveto 3 ph -36 xy rlineto
2 ph 108 xy rlineto stroke
0 0 1 ph -36 0 blue 2 ph 0 0 ph 180 288 red green}
{/n n 1 sub def 0 0 0 kite1 4 ph -36 xy 144 dart1
/n n 1 add def} ifelse
grestore end} def

300 400 translate
/n 3 def init
0 72 288 {/a exch def 0 0 a kite0 0 0 a kite1} for
showpage

これで描いたn=3の図は次のようだ.

ここで, 最初に作る配列 phs は, φの等比級数を保持する. 再帰が深くなる, つまり, nが小さくなるにつれ, phsの添字の小さい方を使うようになっている.

thickとthinも同様なので, 省略する.

余計なことながら, YouTubeにdeflationの面白い動画があった.

2010年8月20日金曜日

都府県の接続グラフ

TAOCP vol.4にContiguous United States of Americaというグラフがある.



アメリカのハワイとアラスカを除く, 本土の48州とDCの隣接関係をグラフにしたものだ. 州名は郵便コードか何かの略号で書いてあるが, 大体は推測がつく. 殆んどは3つの州を頂点とする三角形の辺のグラフで, 注目すべきは左のUT, CD, NM, AZの四角と,中央上のWI, MI, IN, ILの四角だ. 前者は有名な4 sates cornerであり, 後者はミシガン湖である.

前から日本の都府県についても書いてみたいと思っていたところ, やっと書くことが出来た. もちろん北海道と沖縄は入っていない.


関東地方では古河のあたり, 茨城と埼玉が隣接し, 栃木と千葉を隔てているなど, クリティカルな場所がはっきりする. 隣接県の数(次数)が多い県は長野(8)であり, 最小は長崎(1)である. 長野が多いのは 信濃の国は十州に境連ぬる国にして の歌の通りだ.


ついでだが, 隣接する辺は次の通り. 頂点の番号は, 都道府県コードである.

((2 3)(2 5)(3 4)(3 5)(4 5)(4 6)(4 7)(5 6)(6 7)(6 15)(7 8)
(7 9)(7 10)(7 15)(8 9)(8 11)(8 12)(9 10)(9 11)(10 11)(10 15)
(10 20)(11 12)(11 13)(11 19)(11 20)(12 13)(13 14)(13 19)
(14 19)(14 22)(15 16)(15 20)(16 17)(16 20)(16 21)(17 18)
(17 21)(18 21)(18 25)(18 26)(19 20)(19 22)(20 21)(20 22)
(20 23)(21 23)(21 24)(21 25)(22 23)(23 24)(24 25)(24 26)
(24 29)(24 30)(25 26)(26 27)(26 28)(26 29)(27 28)(27 29)
(27 30)(28 31)(28 33)(29 30)(31 32)(31 33)(31 34)(32 34)
(32 35)(33 34)(34 35)(36 37)(36 38)(36 39)(37 38)(38 39)
(40 41)(40 43)(40 44)(41 42)(43 44)(43 45)(43 46)(44 45)
(45 46))

辺の数は86である.

2010年8月10日火曜日

リレー加算器

以下は夏休み工作の記である.

前回の説明通り, 全加算器は4極双投のリレー2個で実現出来る. 田中君によると, 4極双投のリレーは市販されている, というので, 秋葉原へ出かけ, パナソニック電工の12volt駆動のAHJ3241(1個270円)を8個購入した. (写真はパナソニック電工のカタログから)


スイッチは, A0,...,A3, B0,...,B3の入力用8個の他, C0 inも欲しい(最下段に入れる繰上げは, 減算で必要). 始めの8つはオンオフだが, 最後のCは, 切替えなので, 全部同じ切替えスイッチ9個を購入(1個130円). LEDはS0,...,S3の4個の他, C3 outも含めて5個購入(12volt用 1個30円)した. 従って, 合計3480円也であった.




これが組立て用の図である. 下の方に8個のリレーがあり, 上下のリレーが対の全加算器が4個配置されている. 赤が12voltの線. 青がアースである.

上が操作パネルの裏側で, 上の方にLEDが5個. その下にスイッチが二進数A用に4個, B用に4個. 下から入れる繰上げCin用に1個配置する.

スイッチは下向きを0, 上向きを1にしたい. 裏の端子では, スイッチを上向きにすると, 中と下が接続する. 従って, 左に書いてある繰上げ用スイッチは, 上の端子が¬Cin用である.

AのリレーのA0,A1,A2,A3の端子, BのリレーのB0,B1,B2,B3の端子, それぞれ操作パネルの対応するスイッチに, B0リレーの9と10番の端子は, スイッチCinの¬CとCに接続する.

Aの各リレーの10番の端子Sは, 対応するLEDへ, B3リレーの12番の端子Cは, CoutのLEDに接続する.

この図では操作パネルの方が幅狭く描いてあるが, 実際はパネルは幅が100ミリ. リレーは1個の幅が21.5ミリ, 4個で86ミリであり, リレーの方が幅は狭い.

パネルの下の方, 裏面にリレーを固定してある.

電源は, パソコン用のアダプターを利用した.

完成した加算器の写真を以下に示す. (Coutのラベルのところは, LEDの基盤取りつけネジの穴と重なってしまった. LEDの取りつけ方は最後まで決らなかったので, こういうことに.)


表側


裏側

工作は私の属している研究所の工作室で行った. 電気回路をいじるのは実に久しぶりなので, 研究所の宇夫君に最新技術をいろいろ教えて貰った. 工作後一番の感想は, 視力がずいぶん落ちていることだが, 年のことを考えるとまぁ仕方ないか.

2010年8月8日日曜日

十六進乗算表

1985年頃, NTTの電気通信研究所でLispマシンの研究をしていた. そのマシンをELISという. Eは通研 (Electro Communication Laboratories) LISはLispである.

そのプロジェクトで作られた計算機の何台かが, JAIST(北陸先端科学技術大学)にあり, 8月7日から9日にかけて, JAISTでELIS復活祭というイベントが行われた. 20年くらい前の計算機の電源を入れ, もう一度走らせてみようというのである.

同時に多くの発表もあった.

私もLispプログラマなので, お誘いを受け, 参加し, 楽しい時を過している.

そのプロジェクトの最後の時期に TAO/SILENTというのがあり, その発表の時, 当時開発メンバーが十六進法をどう読んでいたかという話題になった. つまりFAB1を「エフエイビーイチ」ではなく, フタブイのように読むのである. これは面白いと思い, 早速(十進の九々に相当する)十六進の乗算表を作ってみた.

まず読み上げ方. 数の次が標準の読み方で, 読みにくい時はかっこ内が使える.

0 レ(オ)
1 イ
2 ニ
3 サ(ミ)
4 ヨ(シ)
5 ゴ
6 ロ(ム)
7 ナ
8 ハ(バヤ)
9 ク
A タ(カ)
B ブ(ボ)
C チャ
D ド
E テ(ケ)
F フ

とりあえず, 十六進数で乗算表を作ると, 次のようになる. 九々も0と1の段がないが, ここでも2の段から始まる. 最初は十六進でもニニガシである. 横に14個並ばないので, 縦長になり過ぎたが.

2204 2306 2408 250A 260C 270E 2810
2912 2A14 2B16 2C18 2D1A 2E1C 2F1E

3206 3309 340C 350F 3612 3715 3818
391B 3A1E 3B21 3C24 3D27 3E2A 3F2D

4208 430C 4410 4514 4618 471C 4820
4924 4A28 4B2C 4C30 4D34 4E38 4F3C

520A 530F 5414 5519 561E 5723 5828
592D 5A32 5B37 5C3C 5D41 5E46 5F4B

620C 6312 6418 651E 6624 672A 6830
6936 6A3C 6B42 6C48 6D4E 6E54 6F5A

720E 7315 741C 7523 762A 7731 7838
793F 7A46 7B4D 7C54 7D5B 7E62 7F69

8210 8318 8420 8528 8630 8738 8840
8948 8A50 8B58 8C60 8D68 8E70 8F78

9212 931B 9424 952D 9636 973F 9848
9951 9A5A 9B63 9C6C 9D75 9E7E 9F87

A214 A31E A428 A532 A63C A746 A850
A95A AA64 AB6E AC78 AD82 AE8C AF96

B216 B321 B42C B537 B642 B74D B858
B963 BA6E BB79 BC84 BD8F BE9A BFA5

C218 C324 C430 C53C C648 C754 C860
C96C CA78 CB84 CC90 CD9C CEA8 CFB4

D21A D327 D434 D541 D64E D75B D868
D975 DA82 DB8F DC9C DDA9 DEB6 DFC3

E21C E32A E438 E546 E654 E762 E870
E97E EA8C EB9A ECA8 EDB6 EEC4 EFD2

F21E F32D F43C F54B F65A F769 F878
F987 FA96 FBA5 FCB4 FDC3 FED2 FFE1

これはこんなプログラムで生成する.

(for-each (lambda (a)
(newline) (newline)
(for-each (lambda (b)
(display (string-append (number->string a 16)
(number->string b 16)
(number->string (quotient (* a b) 16) 16)
(number->string (modulo (* a b) 16) 16) " ")))
(a2b 2 9))
(newline)
(for-each (lambda (b)
(display (string-append (number->string a 16)
(number->string b 16)
(number->string (quotient (* a b) 16) 16)
(number->string (modulo (* a b) 16) 16) " ")))
(a2b 9 16))) (a2b 2 16))

関数 (a2b m n) は, (m m+1 ... n-1) のリストを作る. このプログラムを多少手直しすると, 九々ではなく, フフが得られる. 長い段は途中で折り返してあるので, 要注意.

気づている人もあろうが, 九々では積が1桁の場合, ニニガシのように, 「ガ」を挿入する. 下の表も「ニニガヨ」から始まる.


こんなもの, 使えるだろうか.

ニニガヨ ニサガロ ニヨガハ ニゴガタ ニロガチャ ニナガテ
ニハイレ
ニクイニ ニタイヨ ニブイロ ニチャイハ ニドイタ ニテイチャ
ニフイテ

サニガロ ササガク サヨガチャ サゴガフ サロイニ サナイゴ
サハイハ
サクイブ サタイテ サブニイ サチャニヨ サドニナ サテニタ
サフニド

ヨニガハ ヨサガチャ ヨヨイレ ヨゴイヨ ヨロイハ ヨナイチャ
ヨハニレ
ヨクニヨ ヨタニハ ヨブニチャ ヨチャサレ ヨドサヨ ヨテサハ
ヨフサチャ

ゴニガタ ゴサガフ ゴヨイヨ ゴゴイク ゴロイテ ゴナニサ
ゴハニハ
ゴクニド ゴタサニ ゴブサナ ゴチャサチャ ゴドシイ ゴテシロ
ゴフシブ

ロニガチャ ロサイニ ロヨイハ ロゴイテ ロロニヨ ロナニタ
ロハサレ
ロクサロ ロタサチャ ロブシニ ロチャシハ ロドシテ ロテゴヨ
ロフゴタ

ナニガテ ナサイゴ ナヨイチャ ナゴニサ ナロニタ ナナサイ
ナハサハ
ナクサフ ナタシロ ナブシド ナチャゴヨ ナドゴブ ナテロニ
ナフロク

ハニイレ ハサイハ ハヨニレ ハゴニハ ハロサレ ハナサハ
ハハシレ
ハクシハ ハタゴレ ハブゴハ ハチャロレ ハドロハ ハテナレ
ハフナハ

クニイニ クサイブ クヨニヨ クゴニド クロサロ クナサフ
クハシハ
ククゴイ クタゴタ クブロサ クチャロチャ クドナゴ クテナテ
クフハナ

タニイヨ タサイテ タヨニハ タゴサニ タロサチャ タナシロ
タハゴレ
タクゴタ タタロヨ タブロテ タチャナハ タドハニ タテハチャ
タフクロ

ブニイロ ブサニイ ブヨニチャ ブゴサナ ブロシニ ブナシド
ブハゴハ
ブクロサ ブタロテ ブブナク ブチャハヨ ブドハフ ブテクタ
ブフタゴ

チャニイハ チャサニヨ チャヨサレ チャゴサチャ チャロシハ
チャナゴヨ チャハロレ
チャクロチャ チャタナハ チャブハヨ チャチャクレ チャドクチャ
チャテタハ チャフブヨ

ドニイタ ドサニナ ドヨサヨ ドゴシイ ドロシテ ドナゴブ
ドハロハ
ドクナゴ ドタハニ ドブハフ ドチャクチャ ドドタク ドテブロ
ドフチャサ

テニイチャ テサニタ テヨサハ テゴシロ テロゴヨ テナロニ
テハナレ
テクナテ テタハチャ テブクタ テチャタハ テドブロ テテチャヨ
テフドニ

フニイテ フサニド フヨサチャ フゴシブ フロゴタ フナロク
フハナハ
フクハナ フタクロ フブタゴ フチャブヨ フドチャサ フテドニ
フフテイ

この方のプログラムは, 以下の通り.

(define kana0 '
("レ" "イ" "ニ" "サ" "ヨ" "ゴ" "ロ" "ナ"
"ハ" "ク" "タ" "ブ" "チャ" "ド" "テ" "フ"))
(define kana1 '
("ガ" "イ" "ニ" "サ" "シ" "ゴ" "ロ" "ナ"
"ハ" "ク" "タ" "ブ" "チャ" "ド" "テ" "フ"))
(define (foo n)
(display (list-ref kana1 (quotient n 16)))
(display (list-ref kana0 (modulo n 16))))
(define (bar n)
(display (list-ref kana0 n)))

(for-each (lambda (a)
(newline) (newline)
(for-each (lambda (b)
(bar a) (bar b) (foo (* a b)) (display " "))
(a2b 2 9))
(newline)
(for-each (lambda (b)
(bar a) (bar b) (foo (* a b)) (display " "))
(a2b 9 16)))
(a2b 2 16))

2010年8月2日月曜日

リレー加算器

リレーの話だ. 運動会とは関係ない. ここでいうリレーは電気部品で, 電気信号を継続するという意味から, 和語では継電器という.

私が学生時代に所属していた研究室では, 私が研究室に入る少し前までは, 論理回路の素子はリレーであり, それを使って, ディジタル抵抗計を試作したりしていた. 従って, 研究室にはリレーが沢山あった. 物理教室の実験室には, 直流電源が来ていたから, リレーの実験は簡単で, 私も電話交換機に使う回転スイッチを回して遊んだりしたものだ.

当時, 工学部にはリレー回路の講義があった. 私は大学院生になると, 応用物理学科の磯部教授の, William Keisterの教科書The Design of Switching Circuitsによる講義に出席した. もちろん当時は, その後で私もその学科のメンバーになるとは, まったく思わなかった.

私自身は, ちょうど後藤さんがパラメトロンを発明した直後だったから, 論理回路といえばパラメトロン一点張りであり, リレーの経験はない. そこで一度リレーで二進4桁くらいの加算器を作ってみたいと思っていた.

昨年の夏だったか, 日本橋学館大学の田中君が, リレーの加算器を見せてくれたのに刺激され, 私も実現に向け, 少しずつ進んだ.

まずリレーのおさらいを少々.

リレーには働きにより3種類あり, 左からmake, break, transferという. 下の箱にワイヤが巻いてあるのが電磁石で, 電流が流れると(1の時), その磁力により, 普段(0の時)はバネで上の位置にあるスイッチが下の位置になる.

makeは0の時, 左右の回線が切れ, 1の時, 繋がる. breakは0の時, 回線が繋がり, 1の時, 切れる. transferは, 0の時, 左の回線は右の上と繋がり, 1の時, 下と繋がる. このtransferは, 切替えの時, 左の回線は一旦, 上下のいずれとも切れるが, transferには, 一旦上下両方に繋がるタイプもある.

注意: リレーの回路は, このように, 電流が流れていない時(0の時)の図を描くことになっている.


こういうリレーを想定し, 2ヶ月程前, 二進加算器のリレー回路なんか簡単さと思って, 下のように設計した.


図が左右に長いので, 途中で分割してあり, 上の図の右が, 下の図の左に続く.

一番下の桁a0, b0の和がs0で, 繰上げがc1である. 次のa1, b1, c1の和がs1で, 繰上げがc1である. のように見る.

そして, 田中君に意見を求めたところ, 繰上げにリレーを使うと遅くなるからよくないというコメントが早速来た. そしてZuseの回路(dual-rail-carry full adder)を教えてくれた.



これがニ進1桁分, いわゆる全加算器である. 1桁のaとbと, 下の桁からの繰上cinを足す. この入力に対し, 出力はこの桁の和outと, 上の桁への繰上げcoutである.

この回路の最大の味噌は, 繰上げに裏回路があることだ. (dual circuitという.) つまりcinと¬cin, coutと¬coutがある. 回路図では, notは上線で表現するが, htmlでは書けないので, この説明では¬を使う.

左のb=0と書いてある上が, リレーbで動く(つまりbの値の0,1で動く)リレーで, 4回路のtransferになっている. 4極双投(four pole double throw)ともいう.

右手のa=0と書いてある上がのリレーaである. aとbは上下が反対になっているが, 図を簡単にするためだ.

この他, V+には, 出力を駆動する電源を繋ぐ. (常時1になっている.)

前述のように, aのリレーもbのリレーも0の時の図なので, この時の注目すべき回線を赤線で示す. つまりこの時は, V+の1が¬coutを1にしている. cin, ¬cinはいずれかが1のはずだが, ¬cinが1なら, その先はどこにも繋がっていないので, 関係なし. cinが1なら, 下から繰上げがあり, outを1にする.

全加算器の真理値表は下のようだ. 計算機屋なら夢の中でも書けるはずだ.


入力 出力
a b c s c ¬c
0 0 0 0 0 1
0 0 1 1 0 1
0 1 0 1 0 1
0 1 1 0 1 0
1 0 0 1 0 1
1 0 1 0 1 0
1 1 0 0 1 0
1 1 1 1 1 0

さてZuseの回路のチェックは次のようにする. まずoutが1になる条件を調べる. aとbの組合せを4通り書き, outが1になる条件を探す. つまりoutから逆にV+かcinか¬cinを見るわけである. 0 0ならc=1だ. 0 1ならbの節点を上に切替えてみると, ¬cinに辿り着くからc=0だ. このようにして書くと

a b c
Out= 0 0 1
0 1 0
1 0 0
1 1 1

が得られ, 二進3変数の和であることが分かる. 同様にして, cout, ¬coutも書いてみる. xはそういうcはないこと, oはcに関係なく1になることを示す.

a b c
cout 0 0 x
0 1 1
1 0 1
1 1 o

a b c
cnot 0 0 o
0 1 0
1 0 0
1 1 x

というわけで, この全加算器回路は正しいことが分かる. 繰上げが伝搬する場合, aがすべて0, bがすべて1, 一番下のcinが1なら, 右から来たcinは左のcoutから出ていき, aとbのリレーの動作が終った途端に一番上のcoutを1にする.

ウェブページを探したら, リレーを使って計算機を作ったという, やはり物好きな報告を見つけた.
http://web.cecs.pdx.edu/~harry/Relay/RelayPaper.htm
Konrad Zuseの全加算器の説明の後, "Except for the full adder, I designed all circuits used in my computer."と書いてあるから, この加算器は誰にも真似の出来ない完成品なのであろう.

Processingで描いたこの回路のシミュレータが, http://playground.iijlab.net/~ew/relayadder/relayadder.htmlにある. (ダウンロードには多少時間がかかります.) 右下のa,b,cの枠内をクリックすると, その値の0,1が反転する.
キーボードのキーを押しても値は変わる.

2010年7月29日木曜日

再帰曲線

ドラゴン曲線に関する最初のブログ(2009年6月19日)に書いたことだが,...

左折を+, 右折を-と表示すると
+ + - + + - - + + + - - + - -
になる. 記号は15個あるが, 中央の + が最初の折り目で, その両側はもともと重なっていたのを開いたから, 並び順が中央を中心にして対象で, しかも+と-が入れ替わっている.

これを X2 + Y2 と書くと,
X2 = + + - + + - -
Y2 = + + - - + - -
である. するとこの両枝は同じなので,
X2 = X1 + Y1
Y2 = X1 - Y1
X1 = + + -
Y1 = + - -
最後は
X0 = +
Y0 = -
この+と-のパターンはなんだろうと思っていたが, はたと思いついたのは, Gray codeであった.

下の図は通常の二進法とGray codeの対応を示す. Gray codeの右の+と-は, Gray codeの1の数が, すぐ上のcodeのそれと較べて, 増えたか減ったかを示す. これが上のパターンと同じなのである.




なぜかというと, 枠で囲った部分がX2とY2で, その間が, 一番左の桁が1になったことで + になっているのである.

この図からははみ出すが, 16番のところは右から5ビット目が1になり, +となり, 24番のところは右から4ビット目が0になり, -になるのである.

一方, 通常の二進法の右の小さい丸は, 右から1ビット目と2ビット目, 2ビット目と3ビット目のように, 連続して2つの1が続く状態が現れた場所を示す.

1が連続すると, Gray codeの作り方から分かるように, 排他的論理和をとるので, 1が減るのである.

従って, -が現れるのは, 4n+3か, 8n+6か, 16n+12か, ...である.

これだけ分かるとドラゴン曲線のプログラムは書ける.


(define (test n p)
(cond ((= (modulo n 4) 3) #f)
((= (modulo n 8) 6) #f)
((= (modulo n 16) 12) #f)
((= (modulo n 32) 24) #f)
((= (modulo n 64) 48) #f)
((= (modulo n 128) 96) #f)
((= (modulo n 256) 192) #f)
((= (modulo n 512) 384) #f)
(else #t)))

山場はこのtestだ. こういつまでも書くわけにはいかないから, nがどんなに大きくてもいいように, 下のプログラムでは書き直してある.


(define (moveto x y)
(display (number->string x)) (display " ")
(display (number->string y)) (display " moveto\n"))

(define (rlineto x y)
(display (number->string x)) (display " ")
(display (number->string y)) (display " rlineto\n"))

(define limit 64)

(define (test n p)
(cond ((> p (* 2 n)) #t)
((= (modulo n p) (* 3 (/ p 4))) #f)
(else (test n (* p 2)))))

(define (foo n dx dy)
(rlineto dx dy)
(if (< n limit)
(if (test n 4) (foo (+ n 1) (- dy) dx)
(foo (+ n 1) dy (- dx)))))
(moveto 0 0)
(foo 1 20 0)

2010年7月24日土曜日

再帰曲線

このブログにドラゴン曲線のことを書いたのは, 1年くらい前のことである.

ウェブでドラゴン曲線の画像を眺めていたら, 4つのドラゴン曲線を組合わせて, 平面を埋めている絵があった. 早速描いてみることにした. 次がその戦果である.


もちろんドラゴン曲線は角で直角に曲がるのだが, 通過の仕方が分かるように, 四角の部屋を丸く掃く居候モードで描いてある.

問題は, それぞれの色のドラゴンが重ならないかということだ.

ちゃんと証明するのは面倒だが, ドラゴン曲線をフラクタルで描くステップを考えると, 当たり前のような気がする.

フラクタルでの描き方は次のようだ.



まず0次の線を赤のように引く. 左から右に向って引いている積りだ.
次に1次の線を青のように, 赤の線分を右に三角形に膨らませるように引く.
2次の線は, 橙のように, 青の線分を最初は右に, 次は左に膨らませる.
3次の線は, 緑のように, 橙の線分を右, 左, 右, 左と膨らませる.
これを適当な次数まで繰り返す.

描画アルゴリズムとしては, n次のフラクタルを描きたいというサブルーチンで始める. 次数が0なら始点と終点を引く. そうでないなら, n-1のフラクタルのサブルーチンを呼ぶのである.


このようにして描いた, 0次, 2次, 4次, 6次は次のとおりである. 8次の図はこのブログ先頭のものだ.






重ならない理由がなんとなく納得出来たであろうか.

ついでに10次と12次の図も示すと次のとおり.


2010年7月10日土曜日

ビットの反転と置換

Hacker's DelightにもTAOCPにもSchroeppelのビット反転法の話が登場する. もとはHAKMEMの167項の一部にある話である.

Hacker's Delightでは(102ページ), 64ビットレジスタで7ビットxの反転を((x*0x40100401)&0x442211008)%255 と書いてあり, TAOCP(V4F1,26ページ),0<=x<2gの反転をt←((ax mod 2n)&b, y←((ct)mod 2n)>>(n-g) ただしn=g2, 0≤x<2g, a=(2n+g-1)/(2g+1-1), b=2g-1(2n-1)/(2g-1), c=(2n-g-1)/(2g-1-1) とある.

まずDelightの方.

((x*A)&B)%255
ただし A=0x40100401, B=0x442211008 と書き直す.

このAとBを二進にする.
A=100 0000 0001 0000 0000 0100 0000 0001
つまり4つの1の間に0が9個ある.
B=100 0100 0010 0010 0001 0001 0000 0000 1000

従って 7ビットgfedcbaにAを掛けると

gfedcba
1000000000100000000010000000001 A
gfedcba000gfedcba000gfedcba000gfedcba A*x
10001000010001000010001000000001000 Bでマスクすると
e a f b g c d これらが残る
43210 6543210 6543210 6543210 6543210 8ビット内の位置
11111111 28-1で割った剰余
abcdefg 反転出来た.

ビットに右から番号をつけるとaの位置は0, bは1, ..., gは6.
Aを掛けるとaの位置は0,10,20,30になり, 8の剰余は0,2,4,6になる. Bで6のaを取る. bの位置の8の剰余は1,3,5,7で5のものを取る. そういう仕掛けである.

TAOCPの方は, 定数に(2p-1)/(2q-1)が多いが, これはqビットごとの1のあるパターンを作る常套手段である.


g=3とすると, n=9. 従って

a=100010001
b=100100100
c=10101

x=rqp
これにa=100010001を掛けると
rqp0rqp0rqp
2nでmodをとると
p0rqp0rqp
b=100100100
でマスクする
p00q00r00
c=10101を掛けると
p00q00r00
p00q00r00
p00q00r00
p0pqpqrqr0r00
2nでmodをとると
pqrqr0r00
n-gビット右シフトすると
pqr
と反転出来る.


これは, 反転すべきパターンを位相をずらして複製し, 必要な部分を取り出し, Delightのように2n-1で割った剰余で揃えるか, TAOCPのように(割り算は出来ない)cをかけて混ぜ合わせ, 途中に出来た部分を取り出すかである.


TAOCPのやり方で, もっと一般的な置換が出来ないか考えたのが, 今回のブログのテーマである. 例として64ビットレジスタで8ビットのパターンを任意に置換する.

x=76543210 (ビットの名前)
とする.
置換は (0,1,2,3,4,5,6,7)->(3,2,4,1,6,0,5,7)
つまり75061423にしたい.


次のようにする.
a=0x8040201008040201
b=0xbfdfeff7fbfdfeff
y=(xa % 2^64) & b
c=0x0101010101010101
d=0x4020100804020100
z= ((xc % 2^64) >> 1) & d
m=0x14012000000a4080
(((y|z & m) * c) & (2^64-1)) >> 56

Schemeでは
(define (genperm x m) ;m permutation mask
(let* ((a #x8040201008040201) (b #xbfdfeff7fbfdfeff)
(c #x0101010101010101) (d #x4020100804020100)
(y (band (modulo (* x a) (expt 2 64)) b))
(z (band (>> (* x c) 1) d)))
(>> (band (* (band (bor y z) m) c) (- (expt 2 64) 1)) 56)))
(define m #x14012000000a4080)
実行してみると
(genperm #b11110000 m) => #b11010100
(genperm #b11001100 m) => #b10010011
(genperm #b10101010 m) => #b11001001

yはこうなる.

y= mask
76543210 ff 64ビットの内 最右の8ビット
6543210x fe
543210x7 fd
43210x76 fb
3210x765 f7
210x7654 ef
10x76543 df
0x765432 bf 最左の8ビット x: ドントケアビット

xc>>1 はこうなる.

07654321
07654321
07654321
07654321
07654321
07654321
07654321
?7654321 ?: 右シフトで左から挿入されたビット

z= ((xc mod 2^64) >>1) & d

z= mask
xxxxxxxx 00
xxxxxxx1 01
xxxxxx2x 02
xxxxx3xx 04
xxxx4xxx 08
xxx5xxxx 10
xx6xxxxx 20
x7xxxxxx 40

y|z は各列に0から7を1つずつ含む.

y|z & 0x14012000000a4080 <= 置換用マスク

76543210 7 80
65432101 5 40
54321027 1 2 0a
43210376 00
32104765 00
21057654 0 20
10676543 3 01
07765432 6 4 14
75061423

0x0101010101010101を掛ける.

7
75
75 1 2
75 1 2
75 1 2
750 1 2
75061 23
75061423

modulo 2^64 >> 56 => 75061423

完成!

置換用マスクの作り方.

00 09 18 27 36 45 54 63
08 01 10 19 28 37 46 55
56 17 02 11 20 29 38 47
48 57 26 03 12 21 30 39
40 49 58 35 04 13 22 31
32 41 50 59 44 05 14 23
24 33 42 51 60 53 06 15
16 25 34 43 52 61 62 07

y|zの行列を眺めながら作った上の表は, 0を0へ移動するには0にマスクを置く. 0を1へ移動するには9にマスクを置く. ... 0を7へ移動するには63にマスクを置く. ...7を7へ移動するには7にマスクを置く. のように読む

この表は pをqに移動するマスク位置を計算するプログラムでも計算出来る.

(define (pq p q)
(cond ((< p (+ q 1)) (+ (* -8 p) (* 9 q)))
((= p (+ q 1)) (+ (* 8 p) q))
((> p (+ q 1)) (+ 72 (* -8 p) (* 9 q)))))

上の例 0->5 1->3 2->1 3->0 4->2 5->6 6->4 7->7

(number->string
(apply +
(map (lambda (p q) (expt 2 (pq p q)))
;from pos to pos
'(0 1 2 3 4 5 6 7) '(5 3 1 0 2 6 4 7))
) 16)
=> "14012000000A4080"

2010年7月1日木曜日

再帰曲線

第1回のArtificial LifeのProceedings(1987年)を見ていたら, 下のようなフラクタルの絵があった.



当然どう描くのか不思議であるが, Wikipediaに85度というヒントがあったので, 描いてみたら, 何とも簡単であった.



をまず描くのである. この4本の線をまたフラクタルにすると,





と細かくなり, 何段目かで最初の図になる.



原点からx軸に沿い, 長さdの線を描くとする. その線を途中で折曲げて上のように描きたいから, まずd1を計算する. d1=d/(2+2*cos 85)なのは容易に分かる. 従ってdの線を描く代りに, x軸に沿い, 長さd1の線を描き, 原点を(d1,0)へ移し, そこで85度回転する. また新しいx軸に沿い,長さd1の線を描き, 原点を(d1,0)へ移し, 今度は-170度回転する. 後は図に従い, 同様にやる.

長さdの線を描く代りに, 短い線を何本か描くので, 再帰呼出しになっている. 当然どこかで止めなければならない. nなるパラメータを1つ用意し, 下請けを呼ぶ時, nを1引く. nが0で呼ばれたら, dの直線を直接引く.

PostScriptのプログラムは以下の通り.

/draw {4 dict begin
/n exch def /d exch def %パラメータを取り込む
n 0 gt % n>0なら
{/d1 d 85 cos 1 add 2 mul div def /n1 n 1 sub def
gsave %環境を待避
d1 n1 draw
d1 0 translate
85 rotate
d1 n1 draw
d1 0 translate
-170 rotate
d1 n1 draw
d1 0 translate
85 rotate
d1 n1 draw
grestore} %環境を回復
{0 0 moveto d 0 lineto stroke} ifelse end} def % n=0なら直接描く

40 40 translate
400 6 draw %起動


最終の姿は, 予想もしないものだ.

2010年6月13日日曜日

Ulam数

TAOCP V4F1にUlam数の話題があった(演習問題7.1.3--141).

なんとか数にはFibonacci数, Mersenne数, Goedel数, Meertens数, Erdos数など, 数々あるが, Ulam数はあまりポピュラーではないようで, 私は始めて聞いた. A New Kind of Scienceにはちゃんとあった.

こういう数列である.
<U1, U2,...>=
<1,2,3,4,6,8,11,13,16,18,26,28,36,38,47,48,53,57,62,69,72,77,82,87,97,99,...>
で, U1=1, U2=2, Un(n≥3)は, 0<j<k<n(jとkは違う)のUj+Ukの一通りの和で >Un-1の最小のものという定義である.

U3は, >U2だからまず3が候補で, 3=U1+U2が1通りなので, 3である.

U4は, >U3だからまず4が候補で, 同じU2を使う2+2は許されないが, 4=U1+U3が1通りなので, 4である.

U5は, >U4だからまず5が候補だが, 5は1+4と2+3と2通りで出来るので失格. 6は2+4しかないのでU5=6となる.

U6は, 6の次の7が候補だが, 7は1+6と3+4と2通りで出来るので失格し, 3+5=8になる.

こういう計算をするプログラムをSchemeで正直に書いて実行してみると, Fibonacci数以上に再帰が起き, U10くらいで辛抱出来なくなる.

上の数列は, SICPにあるメモ化に書き直して計算したものである. もちろんメモ化で断然速くなる.

演習問題の趣旨は, これをビット計算で高速化せよというものだ. もともとはBITの論文にあるらしいが, 解答を見ると以下のようになっている.

今回のブログは, このアルゴリズムの謎解きである.

aとbの2つのビット列を使う. nとm=Unまで計算が進んだとき, 0≤s<2mのsについて, aのs番目のビットの意味は, [s=1であるかs=2であるかsは<mの異なるUlam数の1通りの和かである]. bのs番目のビットの意味は, [sは<mの異なるUlam数の2通り以上の和である]である. ただし[述語]は述語が真の時1, 偽の時0とする(Iverson記法).

従ってn=2, m=U2=2の時, a=0111(左から順に, s=0だから上の述語に合わず0; s=1だから1; s=2だから1, s=3はU1+U2の1通りだから1). b=0000(どのsも2通り以上の和にならない)

さてn←n+1の時のアルゴリズムは次のようだ.



実際aとbは次のようになる.

n=3, mはam+1から右へ探して最初の1の添字, つまり0111の1番右のa3から探し始めると, いきなり1だからm=3. これがU(n=)3=3である. Un-1=2, Un=3だから, アルゴリズムの但し書き(where以下)に従い, a4...a5は0, またb4...b5は0だから, a, bは4ビットから6ビットになり, a=011100, b=000000にする. aを0から2までと3から5までに分け, bは3から5を使う.

アルゴリズムでは代入の左辺が(am...a2m-1,bm...b2m-1)のように書いてあるが, これは右辺を計算してからam...a2m-1とbm...b2m-1に代入したいという気持である. つまり

新しいaの3から5は, (100⊕011)&¬000=111.
新しいbの3から5は, (100&011)|000=000.

従ってa=011111, b=000000になる.
次はn=4, m=(a4から右へ探すので)4. 但し書きによりa=01111100, b=00000000になり, aの4から7は (1100⊕0111)&¬0000=1011. bの4から7は(1100&0111)|0000=0100. 従ってa=01111011, b=00000100.

では, これはなにをやったか.

a=011100の時, U1,U2,U3は1通りの和で作れるということであった. aの3から5は100だったが, まずU3の1はこれが1通りの和であることを示す. それにaの0から2の011をxorするのは, 左からの位置sはm+sの和で作れることを示す. 3については3+0はない; 4については3+1がある; 5については3+2があるということだ. すでに1通りで出来るものは, xorでその条件がリセットされてしまう.

一方, bm+sはリセットされる時はasとam+sが共に1なので1になり, つまり2通りで和が作れることを記憶する.

またaに戻り, その後, bs+mの否定とandをとるのは, 一度リセットされたら, あとは別の和で作れることが分かっても, 1に戻さないためである.

またbの方は, 前の値とorをとるから, 一旦2通り以上となれば, ずーっとこのビットは設定されたままになる.

かような次第で, このアルゴリズムを1回使うと, 次のUlam数が1つ得られる.

MIT Schemeにはbit-stringというデータ型があり, and, or, xor, not, stringappend, substringがとれ, またビットを探す関数もあるから, このアルゴリズムを書くのにお誂え向きである. ただしビットの番号は右端が0なので, 上の説明とは添字の向きが反対になる.

3から始め, c個のUlam数を計算するプログラムは以下の通り.

(define a (unsigned-integer->bit-string 4 14))
(define b (unsigned-integer->bit-string 4 0))
(define n 2) (define m 2)

(define (ulam c)
(if (> c 0) (let* ((m1) (z) (a1))
(set! n (+ n 1))
(set! m1 (bit-substring-find-next-set-bit
a (+ m 1) (* m 2)))
(display (list n m1)) (newline) ;nとmのリストを出力する.
(set! z (make-bit-string (* (- m1 m) 2) #f))
(set! a (bit-string-append a z))
(set! b (bit-string-append b z))
(set! m m1)
(set! a1 (bit-string-and
(bit-string-xor (bit-substring a 0 m)
(bit-substring a m (* m 2)))
(bit-string-not (bit-substring b m (* m 2)))))
(bit-substring-move-right! (bit-string-or
(bit-string-and (bit-substring a 0 m)
(bit-substring a m (* m 2)))
(bit-substring b m (* m 2))) 0 m b m)
(bit-substring-move-right! a1 0 m a m)
(ulam (- c 1))) 'ok))

実行結果は

(ulam 10)
(3 3)
(4 4)
(5 6)
(6 8)
(7 11)
(8 13)
(9 16)
(10 18)
(11 26)
(12 28)
=> ok

面白いよね.

2010年6月9日水曜日

楕円反射

前回のブログで, 共焦点二次曲線のことに触れた. 共焦点二次曲線もつい描いてみたくなる優美な図である.

教養学部の学生の時, 最初の製図の課題が二次曲線を描くことであった. この時は曲線群ではなく, 1本の楕円と1本の双曲線だったと思う. 1本とはいえ, 烏口で墨入れして描くのだから大変であった.

楕円を描こうと思うと, 長軸 a と短軸 b を決める, つまり楕円の外側の寸法を押さえることから始めるのが普通であろう. 例えばa = 200, b = 150; アスペクト比4:3の楕円を描くことを考えてみる.

xを-200から適当な増分のステップで200まで増やし, xに対応するyを計算して次のx,yまで線を引く, でまぁ出来るわけだが, x=200, -200の辺りは勾配が大きく, xの増分を小さくしなければならなす, x=0の辺りは勾配は小さく, 増分は大きくしたい. 増分を途中で変えるくらいなら, x=200の辺りは, yを独立変数にした方が楽である.

切り替えはx^2/a^2=1/2, y^2/b^2=1/2の付近でやるのがよいかも知れぬ. かくして1象限分を描くPostScriptのプログラムは

/a 200 def /b 150 def /d 5 def %dは増分
/soly {1 dict begin /x exch def %xからyを解く
1 x x mul a a mul div sub sqrt b mul end} def
/solx {1 dict begin /y exch def %yからxを解く
1 y y mul b b mul div sub sqrt a mul end} def
/rd {d div round d mul} def %増分の倍数で丸める
/x1 {a a mul 2 div sqrt rd} def %切り替え点
/y1 {b b mul 2 div sqrt rd} def %切り替え点

/quadrant {0 b moveto %1象限分を描く
d d x1{/x exch def x x soly lineto} for
y1 d neg 0{/y exch def y solx y lineto} for
stroke} def

これを4回呼び出す.

quadrant 1 -1 scale %第1象限 x軸に対称にする
quadrant -1 1 scale %第4象限 y軸に対称にする
quadrant 1 -1 scale %第3象限 x軸に対称にする
quadrant -1 1 scale %第2象限 y軸に対称にする


そのようにした描いたのが上の図だ. 第1象限の赤い格子は, xを変数で描いた部分と, yを変数で描いた部分を示す. x=140, y=105 辺りで切り替わっている.


しかし, 媒介変数を使った楕円の表現の方が, 一気に書けて嬉しい. 0≤t<2πについて, (a*cos(t),b*sin(t)) の点を次々に結ぶ. 以下の共焦点の二次曲線は媒介変数法で描いてある.

さて, 共焦点の楕円を書くには, まずfを決め, 次に離心率eをパラメータとして決め, それらからaとbを求めることになる.

a=f/e, b=√(a^2-f^2)として, 媒介変数法で描画するわけだ.

双曲線は, やはり媒介変数による表現がある. 楕円は通常のcosとsinであったのに, それがcoshとsinhになる. つまりhyperbolic cosineとhyperbolic sineになるわけで, hyperbolic(双曲線の)という修飾があるのは, これで分かる.

さようなわけで,
sinh x=(e^x-e^(-x))/2,
cosh x=(e^x+e^(-x))/2
を定義しなければならない.

こうして描いた共焦点二次曲線の絵が上の図である. 楕円の方は, eを0.5から0.9まで0.05ステップで変え, 双曲線の方はeを1.1から1.5まで0.05ステップで変えて描いてみた. a=f/eは曲線がx軸と交差する点(中央からの距離)なので, 楕円では, 0.5が円に近い(焦点から遠い)方, 双曲線では, 1.1が放物線に近い(焦点に近い)方である.

2010年6月8日火曜日

楕円反射

楕円の内側が鏡だとして, 鏡の中の, 焦点以外の1点からの光が楕円の壁で反射し, その光がまた壁の別の場所で反射し, それを逐次繰り返すとどうなるか. これは高橋秀俊:「数理と現象」の92ページと113ページに書いてある話だ. そこを読むと, これは森口研でXYプロッタを使って描いてみて, 実験したというので, その図が同書に載っている.

森口研にいた私だが, XYプロッタのグループとは別であった私は, 今頃になって, その絵をどのようにして描くかということに興味を持った. 今ならXYプロッタでなく, PostScriptを使うことになる.

結論からいうと, 思っていたより簡単であった. 今回はそれを報告しよう.



この図は長軸の左端Aからある勾配で右方向へ出発した光がBで楕円に当り, 反射した光が次にCで楕円に当り, ... を10回程度したところである.

楕円上の点(x0, y0)から勾配 n/mで楕円の内側に出発した直線が, 反対側で楕円に突き当たる点の座標(x1,y1)を知りたい.

それには,
x2/a2+y2/b2=1 (楕円)
(x-x0)/m = (y-y0)/n (直線)
を解く.

具体的には, a=175, b=150の楕円で, x0=-175, y0=0からm=3,n=1で出発してみる.

まずWolframAlphaで
solve (x+175)/3=y, x^2/175^2+y^2/150^2=1
とやってみると,
Results:
x=-175, y=0
x=48215/373 = 129.021 y=37800/373 = 101.340
が得られたから, Bの座標(x1,y1)は129.021と101.340である.

私がPostScriptで書いたプログラムでやってみると, 当然xの解は2つ得られるが,x0でない方がx1で, それに対してyを計算するとyも上下2つ得れ, そこで, 直線の式に載っている方を採用する. その結果, これでも129.021439 101.340492 が得られて安心する.

次に反射する角度を得なければならない. まず基本知識として, 焦点から出発した光は, もう1つの焦点へ集まるということから, Bと焦点F0, F1を結ぶ線は, 入射角と反射角が等しいという関係にある. この証明は以下のようにやる. (texで書いたものを張り付けた.)



2本のピンに絡めた輪と鉛筆を使って楕円を描く手法を考えてみても, 上の性質は確からしい.



従ってABF0の角をBF1に足せば, Cへの方向が求まる. つまり
θA=tan-1(y1-y0)/(x1-x0).
θ0=tan-1y1/(x1+f).
θ1=tan-1y1/(x1-f).
θC=(θ0A)+θ1.

それを新しい勾配m,nとし, x1, y1を改めてx0,y0として, いまの手順を繰り返すと下のような図が得られる.



短軸の下の方から, かなり上の方向を狙って出発すると, こういう図も得られる. いずれにしろ, 包絡線は, 楕円か双曲線になるが, それもF0, F1を焦点とする, いわゆる共焦点円錐曲線(Confocal Conics)になる.



このようにして出来た図はきれいに見えるが, 領域を万遍なく埋めて, 最後をぴたりと出発点で止めるには, 試行錯誤が必要であった. 簡単な場合を描いてみると





のようになる. なんかLissajous曲線の楕円版を描いている気分であった. 上の2つのきれいな絵は, 出発角度と繰返し回数を何回も調整して得られたものである. XYプロッタなら, こういう実験に時間が掛ったであろうが, PostScriptでは瞬時に結果が得られ, 楽々に実験を進めることが出来た.

2010年6月1日火曜日

菱形六十面体

WolframAlphaのロゴの菱形六十面体. 最初は立体がよく認識出来なかったが, 模型を作ったり計算したりしているうちに, よく分かるようになった. 目がひとりでに凹凸に反応するようになった.

そうなると, いよいよ計算図学のお出ましで, いつものようにPostScriptで描いて見たくなる.



菱形六十面体を描くには, 上の正十二面体の各面に, 中を窪ませて, 5つの菱形を貼ることになる. その座標を計算しなければならない.

この正十二面体を描いた時の各稜の長さは2だったので, 下の図の0から10までの頂点の座標を, CDの長さを1として求めたい.

CD=1,
AC=tan 54°=1.376,
AB=BD=1/cos 18°=1.051,
AD=1/cos 54°=1.701,
AA'=h=AD tan 31.7175=1.051
までは簡単である.

この面を起して正十二面体に貼るのが次の図である. これはY軸方向から見た図で, XPが正十二面体の1つの面である. PZは長さ1で, 別の方向の稜の半分である. その結果, 1つの面の11個の点の座標が決る. (XPA'の角度atan(1.618/2.618)が実は31.7...°なので, 上の右の図と同じである.)



(define ver
'((2.618 -1 0) ;0
(2 0 0) ;1
(2.618 1 0) ;2
(1.618 1 .618);3
(1.618 1.618 1.618);4
(1 .618 1.618);5
(1 0 2.618) ;6
(1 -.618 1.618);7
(1.618 -1.618 1.618);8
(1.618 -1 .618);9
(1 0 .618) ;10
))
(define fs
'((0 9 10 1) (2 1 10 3) (4 3 10 5) (6 5 10 7) (8 7 10 9)))

verはverticesで頂点の座様である. 0.618が多いが, これは√5=2.236...の小数部の半分である. つまり2.618=(√5+3)/2, 1.618=(√5+1)/2, 0.618=(√5-1)/2である. 0,2,4,6,8の座標は, 最初の正十二面体の0,1,2,3,4の座標と当然合っている.

fsは各面を構成する頂点の番号を一定の方向に回りながら指定する.

一方, 点10の座標をとるのに, 菱形の長い方の対角線が2であり, 点10がzx面上にある, つまりy座標が0なことがわかっているので, xとz座標を連立方程式で解けるはずであり,WolframAlphaで解いてみた.

((3+(sqrt 5))/2-x)^2+1+z^2=(5-(sqrt 5))/2,
(1-x)^2+(((sqrt 5)-1)/2-z)^2=(5-(sqrt 5))/2

これからx=1, z= (√5-1)/ 2 が得られた.

1つの面の座標と面の番号が確定すると, 後はこの座標の鏡像をとったり, 座標軸に対して回転したりすると, 他の面の情報も苦もなく得られる. 2009年2月20日のブログ参照.

これまで正多面体をずいぶん描いたが, それらは凸形であった. 菱形六十面体には窪みがあり, 同じように描けるとは思えない. 今回はZバッファー方同様に, 60ある面を遠方から手前の順に並べ, 遠方から描くようにした. 結構いい加減な発想だが, うまく行かなければまた考え直すことにして, 兎に角やってみる.

最初に出来た図は下のようである. 計算でこういう絵が出来ると, やはり感動する.

この図は最初の正十二面体と同じ方向から見たように描いてある. 赤線のように飛び出した頂点を結んでみると, 最初の正十二面体に内接していることが理解出来る.

前々回のブログにあったように面に色を塗るとこのようになる.

面白いのは右下の2つの小さい緑の面で, 上の正十二面体では, 右下の破線で囲まれているむこう向きの面である. 手前の2本の実線の稜線が引っ込んだことと, むこう向きの面の中央が窪んだため, 奥の面がこちら向きになって, 見えるようになった.

正面上方の赤い面の方向から眺めた絵を描くと, WolframAlphaの(色違いの)ロゴが出来る. よく見るとロゴの中の星は下が尖っているから, 上下逆さまのようだが, まぁよかろう.

しかし, このロゴを描くだけなら, 遥かに簡単だ. 線の長さはみな同じ. 線の方向も5通りしかない. 線同士の接続順さえ押さえれば, このロゴは描くことが出来る.

/un 80 def %unit length
/ev {36 mul 72 sub dup sin un mul exch cos un mul} def
/rl {ev rlineto} def /rm {ev rmoveto} def
/a {0 rl} def /b {1 rl} def /c {2 rl} def /d {3 rl} def
/e {4 rl} def /f {5 rl} def /g {6 rl} def /h {7 rl} def
/i {8 rl} def /j {9 rl} def /l {1 rm} def /n {3 rm} def
/p {5 rm} def /r {7 rm} def /t {9 rm} def
/u {1 0 0 setrgbcolor} def /v {0 0 1 setrgbcolor} def
/w {0 1 0 setrgbcolor} def /x {1 0.5 0 setrgbcolor} def
/y {0 setgray} def /z {0 0 moveto} def
300 300 translate
/as {z l n g f i h a j c b e d} def %center
u as fill y as z b z d z f z h z j stroke
/k {z r f g a i b j d e} def %6 o'clock
w k fill y k z r g l i stroke
/m {z p d e i g j h b c} def %8 o'clock
x m fill y m z p e t g stroke
/o {z n b c g e h f j a} def %11 o'clock
v o fill y o z n c r e stroke
/q {z l j a e c f d h i} def %1 o'clock
x q fill y q z l a p c stroke
/s {z t h i c a d b f g} def %4 o'clock
v s fill y s z t i n a stroke

この程度に単純に描ける. 周囲の5つのパターンは, 72°ずつ回転しても描けるが, a→c→e→g→i→a, b→d→f→h→j→b, l→n→p→r→t→lの置換が面白く, このままにしてある.

なお, Processingで描いた動画がhttp://playground.iijlab.net/~ew/rhomb/rhomb.htmlにある. マウスを押すと初期状態に戻り, どれかキーを押すと停止する. 理由は不明だが, ロードに1分少々時間がかかるので, 辛抱して待って欲しい.