2010年5月17日月曜日

菱形六十面体

MathWorldのRhombic Hexecontahedronの項に気になる記述があった.

20個の菱形多面体(golden rhombohedra)を合わせると菱形六十面体が作れるというのである. 確かに外から見ると, それぞれの突起は菱形多面体のようである. しかも正十二面体の頂点は20個だから, 数的にも丁度合う.

問題は菱形多面体のとがった頂点の立体角で, 20個合わせて4 πになるかである. 黄金菱形のとがった角の半分はθ=atan(1/φ) (φは黄金比 (√5+1)/2)で, これを3枚合わせた立体角を計算したい. 球面三角法の図を使おう.



この図はOが菱形多面体のとがった頂点で, BOAが黄金菱形のとがった角, Cは底面の黄金菱形の半分で, ∠C=直角になっている. 小文字はOの周りの中心角, 大文字は2面角(dihedral angle)である.

この他 分かっているのは, b=θ; c=2θ; A=2Bの関係である.
sin C/sin c=sin B/sin bから, Bが分かり, Aが分かる.

またMathWorldのSpherical Triangleの項によれば, 2面角A,B,Cで囲まれる球面上の面積はR2[(A+B+C)-π]なので, これで立体角が得られる.

計算すると(単位はradian)

(define phi (/ (+ (sqrt 5) 1) 2)) => 1.618033988749895
(define b (atan (/ 1 phi))) => .5535743588970452
(=31.717474411461005°)
(define c (* 2 b)) => 1.1071487177940904
(define pi (* 4 (atan 1))) =>3.141592653589793
(define C (/ pi 2))
(define B (asin (/ (* (sin C) (sin b)) (sin c))))
=>.6283185307179586

あ, ∠Bはπ/5(=36°)なのか. 2Bの角が5個集まって360°になるわけだから, そういう考えから出発すれば, 何も計算することはなかったのか:-<.

従って2面角Aは2π/5; 球面上の面積は 6π/5-π=π/5.

たしかにこの立体角を20個合わせれば4πになるわ.

2010年5月16日日曜日

菱形六十面体

Wolfram Alphaというすごいサーバーがある. これを楽しんでいるうちに, ひと昔前のエイズウィルスの漫画でよく見たような真っ赤なロゴが気になり, 調べたところ, rhombic hexecontahedronであるということが分かった. rhombicは「菱形の」であり, hexecontaは「六十」なので, 日本語では「菱形六十面体」である.

菱形といっても, WolframMathWorldによれば, golden rhombi(rhombiはrhombusの複数形)という, そんじょそこらの菱形ではない.

つまり長短の対角線が黄金比になっているのである. 菱形六十面体は, 正十二面体の1つの正5角形を5つの黄金菱形で構成した凹面で置換えたものであった.


Wolfram Alphaのページ
には展開図があった. そこで早速PostScriptで写し取り, 印刷した紙を切り抜いて多面体を作る. 菱形内の番号は作業用で, 菱形を次々繋げて描いた順が見て取れる.



PostScriptでの描き方だが, 最後に描いた菱形に対し, 上下左右のいづれかの辺の鏡像を作るよう, 4つのサブルーチンを用意し, それを次々に呼び出している. 大きい枝分かれのあるところは, gsaveとgrestoreを使ってバックトラックする.

この図から作った模型が下の写真である.

こういう模型は最後の糊付けが最大の問題で, この模型も実は最後は留めていない.

一旦出来てみると, 意外に飛び離れた菱形が同一5角形のメンバーなことが分かる. 後のために展開図を塗り分けたが, 4色ではどの面に集まるかの情報が不足なので, 正十二面体の展開図にも番号を入れた.




下の図は, 正十二面体1つの面の構成する凹面5角形を考えてみたものである.

角DB'Eは黄金菱形なので, 黄金比φ=1.61に対して2 cot-1(φ)=63.435°である. その一端を持ち上げると, 角DBEは少しずつ広がり, 72°になるわけだ. その時の傾き角を計算すると, 驚いたことに31.7175°であった. つまり側面図のABと平面図のAD'は平行だった. こうして得られたBCは, 中央の窪みからの頂点の高さである. またD,Eの高さはBCの半分である.

そうと分かれば, 凹面5角形が描ける. 次の図は, いろいろな方向から凹面5角形を描いたもので, 左上が真横から, 続いて視点を15°ずつ増やし, 最後の右下は真上から見たものである.



これと正十二面体のプログラムを合わせれば, 菱形六十面体の絵も描けると思うがそれはいずれまたということにしよう.

2010年5月13日木曜日

オランダ国旗の問題

オランダのEindhoven工科大学にいたDijkstraは, 次のようなプログラミングの問題を出したといわれる.

配列a[0]からa[n-1]に青, 白, 赤いずれかがランダムにある. 1回のスキャンで青を0側, 赤をn-1側, 白を中央に揃えよ.

解法

ポインタb, x, rをb←0, x←n, r←n; と初期化し,
while b≠x do
if a[x-1]=青 then swap(a[b],a[x-1]); b←b+1;
if a[x-1]=白 then x←x-1;
if a[x-1]=赤 then swap(a[r-1],a[x-1]); x←x-1; r←r-1;
のようにループを回す.

次の図から明らかであろう.



0≤ <bの添字は青, x≤ <rの添字は白, r≤ <nの添字は赤, そしてb≤ <xの間は未処理 になっているのがこのループのinvariantである.

この問題(Dutch National Flag Problem)は, 同じ3色の国旗を持つフランス人が, オランダ人のDijkstraにプレゼンとしたといわれるが, そのフランス人の名前は失念した.

2010年5月10日月曜日

Life Game

前回のブログの最後. Schemeの実装はTAOCPのアルゴリズムそのもののコピーで, 実はもっと関数プログラムらしく書きたかった. 書き直したのが以下のプログラムである. ある変数を後段で1回しか使わぬなら, 直接使う場所の書き込んである.

(define (b x- x x+)
(define (and a b) (fix:and a b))
(define (or a b) (fix:or a b))
(define (xor a b) (fix:xor a b))
(define (lsh a b) (fix:lsh a b))
(let* ((a0 (and x- x+)) (b0 (xor x- x+)) (c0 (xor x b0))
(d0 (lsh c0 -1)) (c1 (lsh c0 1)) (e0 (xor c1 d0))
(f1 (or (and b0 e0) (and c1 d0))) (c4 (or (and x b0) a0))
(b1 (lsh c4 1)) (c5 (lsh c4 -1)))
(and (xor (or (and b1 c5) (or a0 f1))
(or (and a0 f1) (or b1 c5)))
(or (xor b0 e0) x))))

このプログラムは, fixnumでも動くが, fixnumは24ビットしかないので, 実用的には長さが自由なbit-string型の方を使いたい. それに対処すべく, 上のプログラムでは, and, or等は別に定義してある.

MIT Schemeのreference manualでBit Stringを見ると, bit-string-and, bit-string-or, bit-string-xorはあるが, シフトはない. そこでシフトはbit-substringとbit-string-append使って実装することになる.

(define (bit-string-lsh a c)
(if (>= c 0)
(bit-string-append (make-bit-string c #f)
(bit-substring a 0 (- (bit-string-length a) c)))
(bit-string-append
(bit-substring a (- c) (bit-string-length a))
(make-bit-string (- c) #f))))

(define bs (unsigned-integer->bit-string 12 #b111010010111))
=>#*111010010111 ;12ビットのテスト用bit stringを作る.
(bit-string-lsh bs 2) => #*101001011100 ;左シフト
(bit-string-lsh bs -2) => #*001110100101 ;右シフト

こんな具合いである.

2010年5月9日日曜日

Life Game

前回の算法はセル1個ごとに次ステップの値を計算するものであった.

通常の計算機は, 1ビットCPUのconnection machineと違って, 語の処理が出来るから, セルの1行について計算したいとだれでも思うであろう.

前回の演習問題の後半はその行ごと処理である.

ある行xと上の行x-と下の行x+からその行の次のステップの値が計算したい. 出来そうでもあり, 無理そうでもある. アルゴリズムは次の通り.

a←x-&x+ (=z3), b←x-⊕x+ (=z4), c←x⊕b, d←c»1 (=z6),
c←c«1 (=z2), e←c⊕d, c←c&d, f←b&e, f←f | c (=z7), e←b⊕e (=z8), c←x&b, c←c | a, b←c«1 (=z5), c←c»1 (=z1), d←b&c, c←b | c,
b←a&f, f←a | f, f←d | f, c←b | c,
f←f⊕c (=S1(z1,z3,z5,z7)), e←e | x, f←f&e.
これをグライダーでやってみよう.


図のように, このグライダーは右下へ進む.
  001000 x-
  000100 x
  011100 x+
  000000
xの行の計算は次のように進む.
a 001000 x-&x+ (=z3)
b 010100 x-⊕x+ (=z4)
c 010000 x⊕b (=xnw⊕xw⊕xsw)
d 001000 c»1 (=z6)
c 100000 c«1 (=z2)
e 101000 c⊕d (=z2⊕z6)
c 000000 c&d (=z2&z6)
f 000000 b&e (=z4&(z2⊕z6))
f 000000 f | c ((z4&(z2⊕z6))∨(z2&z6)=z7)
e 111100 b⊕e (=z4⊕z2⊕z6=z8)
c 000100 x&b (=x&(x-⊕x+))
c 001100 c | a (=(x&(x-⊕x+))∨(x-&x+))
b 011000 c«1 (=z5)
c 000110 c»1 (=z1)
d 000000 b&c (S1の計算開始 b=z5 c=z1 a=z3 f=z7)
c 011110 b | c
b 000000 a&f
f 001000 a | f
f 001000 d | f
c 011110 b | c
f 010110 f⊕c (=S1(z1,z3,z5,z7)
e 111100 e | x (= x∨ z8)
f 010100 f&e

1ステップ後の配列は以下の通り.
  000000
  010100
  001100
  001000

MIT Schemeのfixnum演算を使った実装は以下の通り.

(define (b x- x x+)
(let* ((a) (b) (c) (d) (e) (f))
(set! a (fix:and x- x+)) (set! b (fix:xor x- x+))
(set! c (fix:xor x b)) (set! d (fix:lsh c -1))
(set! c (fix:lsh c 1)) (set! e (fix:xor c d))
(set! c (fix:and c d)) (set! f (fix:and b e))
(set! f (fix:or f c)) (set! e (fix:xor b e))
(set! c (fix:and x b)) (set! c (fix:or c a))
(set! b (fix:lsh c 1)) (set! c (fix:lsh c -1))
(set! d (fix:and b c)) (set! c (fix:or b c))
(set! b (fix:and a f)) (set! f (fix:or a f))
(set! f (fix:or d f)) (set! c (fix:or b c))
(set! f (fix:xor f c)) (set! e (fix:or e x))
(set! f (fix:and f e)) f))

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))))