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分少々時間がかかるので, 辛抱して待って欲しい.