Ulam数

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

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

こういう数列である.
<U₁, U₂,...>=
<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,...>
で, U₁=1, U₂=2, U_n(n≥3)は, 0<j<k<n(jとkは違う)のU_j+U_kの一通りの和で >U_n-1の最小のものという定義である.

U₃は, >U₂だからまず3が候補で, 3=U₁+U₂が1通りなので, 3である.

U₄は, >U₃だからまず4が候補で, 同じU₂を使う2+2は許されないが, 4=U₁+U₃が1通りなので, 4である.

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

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

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

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

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

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

aとbの2つのビット列を使う. nとm=U_nまで計算が進んだとき, 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=U₂=2の時, a=0111(左から順に, s=0だから上の述語に合わず0; s=1だから1; s=2だから1, s=3はU₁+U₂の1通りだから1). b=0000(どのsも2通り以上の和にならない)

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



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

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

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

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

従ってa=011111, b=000000になる.
次はn=4, m=(a₄から右へ探すので)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の時, U₁,U₂,U₃は1通りの和で作れるということであった. aの3から5は100だったが, まずU₃の1はこれが1通りの和であることを示す. それにaの0から2の011をxorするのは, 左からの位置sはm+sの和で作れることを示す. 3については3+0はない; 4については3+1がある; 5については3+2があるということだ. すでに1通りで出来るものは, xorでその条件がリセットされてしまう.

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

またaに戻り, その後, b_s+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

面白いよね.

楕円反射

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

教養学部の学生の時, 最初の製図の課題が二次曲線を描くことであった. この時は曲線群ではなく, 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が放物線に近い(焦点に近い)方である.

楕円反射

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

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

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



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

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

それには,
x²/a²+y²/b²=1 (楕円)
(x-x₀)/m = (y-y₀)/n (直線)
を解く.

具体的には, a=175, b=150の楕円で, x₀=-175, y₀=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の座標(x₁,y₁)は129.021と101.340である.

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

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



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



従ってABF₀の角をBF₁に足せば, Cへの方向が求まる. つまり
θ_A=tan^-1(y₁-y₀)/(x₁-x₀).
θ₀=tan^-1y₁/(x₁+f).
θ₁=tan^-1y₁/(x₁-f).
θ_C=(θ₀-θ_A)+θ₁.

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



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



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





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

菱形六十面体

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

菱形六十面体

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で囲まれる球面上の面積はR²[(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πになるわ.

菱形六十面体

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°ずつ増やし, 最後の右下は真上から見たものである.



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

オランダ国旗の問題

オランダの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にプレゼンとしたといわれるが, そのフランス人の名前は失念した.