再帰曲線

東京大学新聞に2012年度後期日程入試問題が掲載されていた. その総合科目IIにフラクタルの絵があった. 問題はフラクタル次元を計算するものだが, 私としてはその絵を描くプログラムに関心があった.

まず図はこういうものだ.


左からF₁, F₂, F₃

左からG₁, G₂, G₃
問題に記述はこうだ.

はじめに, 1辺の長さ1の正方形をF₀とする. F₀を, 1辺の長さ1/3の9個の小正方形に分割し, 中央の小正方形を取り去ったものをF₁とする.  また, それぞれの小正方形を1辺の長さ1/3のユニットと呼ぶ. 次に, F₁ を構成する8個のユニットのそれぞれについて, これを9個の小正方形に分割し, そこから中央の小正方形を取り去ったものをF₂とする.  また, この分割で得られたそれぞれの小正方形を1辺の長さ1/9のユニットと呼ぶ. 以下, 同様の操作を繰り返して得られる図形をF₃, F₄, ...とする.

PostScriptのこのプログラムは次のようだ.

/n 5 def
/l 540 def
/xs [0 1 2 0 2 0 1 2] def
/ys [0 0 0 1 1 2 2 2] def
/a {n 0 eq {0 0 moveto l 0 rlineto 0 l rlineto l neg 0 rlineto
 closepath fill}
 {/n n 1 sub def
 gsave 1 3 div dup scale
 0 1 7 {/i exch def
 xs i get l mul ys i get l mul gsave translate a grestore} for
 grestore
 /n n 1 add def} ifelse} def
 40 40 translate a

まず全体の次数nを決める. この例では5. 1辺の長さlも決める. 次に8個の正方形を描くので, その順に左下のxとy座標のリストxs, ysを書く. そしてユニットを描く手続きaの定義だ.



n=0なら, 長さlの正方形を書く. そうでないならnを1減らし, スケールを1/3にする. dupはx座標とy座標の両方を1/3にするためのもの. gsaveとそれに見合うgrestoreは, スケールや原点を変えるとき, 以前の環境をスタックするものである. xsとysから各小正方形の原点の座標をとり, 図を各手続きaを呼ぶ.

こうしてn=5を描いたのがこれだ.



次のは少し手ごわい.

はじめに, 1辺の長さ1の正方形G₀を, 1辺の長さ1/2の正方形ユニット4個に分け, 右上のユニットについては, 左下の1辺1/4の正方形を取り去る. こうして得られた図形をG₁とする. 次に, 残った3個の, 欠陥のないユニットのそれぞれについて, 同様の操作を行なう. すなわち, それぞれを4個の1辺の長さ1/4のユニットに分け, 右上にユニットについては左下の1辺の長さ1/8の正方形を取り去る. 得られた図形G₂とする. 続いて, G₂に含まれる1辺の長さ1/4の, 欠陥のない各正方形ユニットについて, 同様の操作を行なう. 得られた図形をG₃とする. これを繰り返して得られる図形の列をG₁, G₂, G₃, ...とする.

目で見ると簡単だが, 欠陥のない正方形をどう判定するか.

私の考え方はこうだ.

各正方形を4区画に分け, 左下, 右下, 左上, 右上の順に下請けに渡すとする. G₀を描く手続きをaとすると, 1/2のサイズで, a,a,aと3回呼び, 最後に左下を取り除いた正方形を描く手続きbを呼ぶ.

aはnを1減らして, a,a,a,bと呼べば良い. bはnを1減らし, 左下を除いて, a,a,aと呼べば良い. そう考えて書いたのが, このプログラムだ. 上のプログラムが理解出来れば, こちらも分かるだろう.

/n 7 def
/l 540 def /l2 l 2 div def
/a {n 0 eq{0 0 moveto l 0 rlineto 0 l rlineto l neg 0 rlineto
 closepath fill}
 {/n n 1 sub def
 gsave
 0.5 dup scale
 gsave 0 0 translate a grestore
 gsave l 0 translate a grestore
 gsave 0 l translate a grestore
 gsave l l translate b grestore
 grestore
 /n n 1 add def} ifelse} def
/b {n 0 eq{
 l2 0 moveto l2 0 rlineto 0 l rlineto l neg 0 rlineto
 0 l2 neg rlineto l2 0 rlineto closepath fill}
 {/n n 1 sub def
 gsave 0.5 dup scale
 gsave l 0 translate a grestore
 gsave 0 l translate a grestore
 gsave l l translate a grestore
 grestore
 /n n 1 add def} ifelse} def
 40 40 translate a

この結果, n=7で描いたのが, 下の図である.



PostScripでこういう図を描くのは思ったより易しい.

再帰曲線

3次元のHilbert曲線. 前回で2次までは描けることが分かった. まぁしかし, 3次を描くとどうなりそうかにも興味があってやってみた. ご用とお急ぎの方のために, 先ず結果をお目にかけるとこのようになる.



中央少し下の赤い点が出発点, 一番上の緑の点が到着点である. 今回はこの絵の描き方を説明したい.

天下り的だが,

(define (z+y^ d x y z)
 (x+z^ (/ d 2) x y z) (y+z^ (/ d 2) (+ x d) y z)
 (y+z^ (/ d 2) (+ x d) (+ y d) z)
 (z+xv (/ d 2) x (+ y d) z) (z+xv (/ d 2) x (+ y d) (+ z d))
 (y-zv (/ d 2) (+ x d) (+ y d) (+ z d))
 (y-zv (/ d 2) (+ x d) y (+ z d)) (x-zv (/ d 2) x y (+ z d)))

は1辺の長さdのz+y^をx,y,zを起点にして描く起動関数である. 描くといっても, 下請け関数が節点の座標を出力するだけで, 実際に描くのはPostScriptにまかせる.

この関数は下請けを8回呼ぶ. 最初のx+z^の場合は, 下請けの起点は貰ってきたx,y,zである. 次のy+z^の場合は, xを(+ x d), つまりxの位置をdだけずらす. z+y^のGrayコードが000, 001, 011,...と始まったのを反映している. だからその次はyも(+ y d)として呼ぶ. このようにして8回呼ぶと終わる.

その下請けは, やはりそれぞれのGaryコードに基づき,

(define (x+z^ d x y z)
 (p x y z) (p x (+ y d) z) (p x (+ y d) (+ z d)) 
 (p x y (+ z d)) (p (+ x d) y (+ z d))
 (p (+ x d) (+ y d) (+ z d)) (p (+ x d) (+ y d) z)
 (p (+ x d) y z))

(define (y+z^ d x y z)
 (p x y z) (p (+ x d) y z) (p (+ x d) y (+ z d))
 (p x y (+ z d)) (p x (+ y d) (+ z d))
 (p (+ x d) (+ y d) (+ z d))
 (p (+ x d) (+ y d) z) (p x (+ y d) z))
...

のようにして, x+z^, y+z^, z+xv, y-zv, x-zvの5個用意する. (p x y z)は(display (list x y z))である. 座標の8個出力後の改行もpの担当である. (z+y^ 2 0 0 0)と起動関数を呼ぶと, 64個の座標

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

が得られ, この座標に適当な倍率を掛けて, PostScriptで描くと


この絵は2次だが, さて, 3次まで描くには, 上にあった起動の関数と下請け関数を統合し,

x+y^ x+yv x+z^ x+zv x-y^ x-yv x-z^ x-zv
y+z^ y+zv y+x^ y+xv y-z^ y-zv y-x^ y-xv
z+x^ z+xv z+y^ z+yv z-x^ z-xv z-y^ z-yv

の24個について関数を書く必要がある. もちろん1個ずつ手書きするわけにはいかない. 当然関数発生関数を書くことになる. それを使い,

(3dgen 2 0 1 0) =>
(define (z+y^ d x y z) (if (= d 1) (begin (p x y z)
(p (+ x d) y z) (p (+ x d) (+ y d) z) (p x (+ y d) z)
(p x (+ y d) (+ z d)) (p (+ x d) (+ y d) (+ z d))
(p (+ x d) y (+ z d)) (p x y (+ z d)))
(begin (x+z^ (/ d 2) x y z) (y+z^ (/ d 2) (+ x d) y z)
(y+z^ (/ d 2) (+ x d) (+ y d) z) (z+xv (/ d 2) x (+ y d) z)
(z+xv (/ d 2) x (+ y d) (+ z d))
(y-zv (/ d 2) (+ x d) (+ y d) (+ z d))
(y-zv (/ d 2) (+ x d) y (+ z d)) (x-zv (/ d 2) x y (+ z d)))))

のような関数を24個揃える. そして(z+y^ 4 0 0 0)と起動すると, 512個の座標

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

(4 0 0)(4 0 1)(4 1 1)(4 1 0)(5 1 0)(5 1 1)(5 0 1)(5 0 0)
(6 0 0)(7 0 0)(7 1 0)(6 1 0)(6 1 1)(7 1 1)(7 0 1)(6 0 1)
(6 0 2)(7 0 2)(7 1 2)(6 1 2)(6 1 3)(7 1 3)(7 0 3)(6 0 3)
(5 0 3)(5 0 2)(4 0 2)(4 0 3)(4 1 3)(4 1 2)(5 1 2)(5 1 3)
(5 2 3)(5 2 2)(4 2 2)(4 2 3)(4 3 3)(4 3 2)(5 3 2)(5 3 3)
(6 3 3)(7 3 3)(7 2 3)(6 2 3)(6 2 2)(7 2 2)(7 3 2)(6 3 2)
(6 3 1)(7 3 1)(7 2 1)(6 2 1)(6 2 0)(7 2 0)(7 3 0)(6 3 0)
(5 3 0)(5 3 1)(5 2 1)(5 2 0)(4 2 0)(4 2 1)(4 3 1)(4 3 0)

(4 4 0)(4 4 1)(4 5 1)(4 5 0)(5 5 0)(5 5 1)(5 4 1)(5 4 0)
(6 4 0)(7 4 0)(7 5 0)(6 5 0)(6 5 1)(7 5 1)(7 4 1)(6 4 1)
(6 4 2)(7 4 2)(7 5 2)(6 5 2)(6 5 3)(7 5 3)(7 4 3)(6 4 3)
(5 4 3)(5 4 2)(4 4 2)(4 4 3)(4 5 3)(4 5 2)(5 5 2)(5 5 3)
(5 6 3)(5 6 2)(4 6 2)(4 6 3)(4 7 3)(4 7 2)(5 7 2)(5 7 3)
(6 7 3)(7 7 3)(7 6 3)(6 6 3)(6 6 2)(7 6 2)(7 7 2)(6 7 2)
(6 7 1)(7 7 1)(7 6 1)(6 6 1)(6 6 0)(7 6 0)(7 7 0)(6 7 0)
(5 7 0)(5 7 1)(5 6 1)(5 6 0)(4 6 0)(4 6 1)(4 7 1)(4 7 0)

(3 7 0)(2 7 0)(2 7 1)(3 7 1)(3 6 1)(2 6 1)(2 6 0)(3 6 0)
(3 5 0)(3 4 0)(3 4 1)(3 5 1)(2 5 1)(2 4 1)(2 4 0)(2 5 0)
(1 5 0)(1 4 0)(1 4 1)(1 5 1)(0 5 1)(0 4 1)(0 4 0)(0 5 0)
(0 6 0)(1 6 0)(1 7 0)(0 7 0)(0 7 1)(1 7 1)(1 6 1)(0 6 1)
(0 6 2)(1 6 2)(1 7 2)(0 7 2)(0 7 3)(1 7 3)(1 6 3)(0 6 3)
(0 5 3)(0 4 3)(0 4 2)(0 5 2)(1 5 2)(1 4 2)(1 4 3)(1 5 3)
(2 5 3)(2 4 3)(2 4 2)(2 5 2)(3 5 2)(3 4 2)(3 4 3)(3 5 3)
(3 6 3)(2 6 3)(2 6 2)(3 6 2)(3 7 2)(2 7 2)(2 7 3)(3 7 3)

(3 7 4)(2 7 4)(2 7 5)(3 7 5)(3 6 5)(2 6 5)(2 6 4)(3 6 4)
(3 5 4)(3 4 4)(3 4 5)(3 5 5)(2 5 5)(2 4 5)(2 4 4)(2 5 4)
(1 5 4)(1 4 4)(1 4 5)(1 5 5)(0 5 5)(0 4 5)(0 4 4)(0 5 4)
(0 6 4)(1 6 4)(1 7 4)(0 7 4)(0 7 5)(1 7 5)(1 6 5)(0 6 5)
(0 6 6)(1 6 6)(1 7 6)(0 7 6)(0 7 7)(1 7 7)(1 6 7)(0 6 7)
(0 5 7)(0 4 7)(0 4 6)(0 5 6)(1 5 6)(1 4 6)(1 4 7)(1 5 7)
(2 5 7)(2 4 7)(2 4 6)(2 5 6)(3 5 6)(3 4 6)(3 4 7)(3 5 7)
(3 6 7)(2 6 7)(2 6 6)(3 6 6)(3 7 6)(2 7 6)(2 7 7)(3 7 7)

(4 7 7)(4 7 6)(4 6 6)(4 6 7)(5 6 7)(5 6 6)(5 7 6)(5 7 7)
(6 7 7)(7 7 7)(7 6 7)(6 6 7)(6 6 6)(7 6 6)(7 7 6)(6 7 6)
(6 7 5)(7 7 5)(7 6 5)(6 6 5)(6 6 4)(7 6 4)(7 7 4)(6 7 4)
(5 7 4)(5 7 5)(4 7 5)(4 7 4)(4 6 4)(4 6 5)(5 6 5)(5 6 4)
(5 5 4)(5 5 5)(4 5 5)(4 5 4)(4 4 4)(4 4 5)(5 4 5)(5 4 4)
(6 4 4)(7 4 4)(7 5 4)(6 5 4)(6 5 5)(7 5 5)(7 4 5)(6 4 5)
(6 4 6)(7 4 6)(7 5 6)(6 5 6)(6 5 7)(7 5 7)(7 4 7)(6 4 7)
(5 4 7)(5 4 6)(5 5 6)(5 5 7)(4 5 7)(4 5 6)(4 4 6)(4 4 7)

(4 3 7)(4 3 6)(4 2 6)(4 2 7)(5 2 7)(5 2 6)(5 3 6)(5 3 7)
(6 3 7)(7 3 7)(7 2 7)(6 2 7)(6 2 6)(7 2 6)(7 3 6)(6 3 6)
(6 3 5)(7 3 5)(7 2 5)(6 2 5)(6 2 4)(7 2 4)(7 3 4)(6 3 4)
(5 3 4)(5 3 5)(4 3 5)(4 3 4)(4 2 4)(4 2 5)(5 2 5)(5 2 4)
(5 1 4)(5 1 5)(4 1 5)(4 1 4)(4 0 4)(4 0 5)(5 0 5)(5 0 4)
(6 0 4)(7 0 4)(7 1 4)(6 1 4)(6 1 5)(7 1 5)(7 0 5)(6 0 5)
(6 0 6)(7 0 6)(7 1 6)(6 1 6)(6 1 7)(7 1 7)(7 0 7)(6 0 7)
(5 0 7)(5 0 6)(5 1 6)(5 1 7)(4 1 7)(4 1 6)(4 0 6)(4 0 7)

(3 0 7)(3 0 6)(2 0 6)(2 0 7)(2 1 7)(2 1 6)(3 1 6)(3 1 7)
(3 2 7)(3 3 7)(2 3 7)(2 2 7)(2 2 6)(2 3 6)(3 3 6)(3 2 6)
(3 2 5)(3 3 5)(2 3 5)(2 2 5)(2 2 4)(2 3 4)(3 3 4)(3 2 4)
(3 1 4)(3 1 5)(3 0 5)(3 0 4)(2 0 4)(2 0 5)(2 1 5)(2 1 4)
(1 1 4)(1 1 5)(1 0 5)(1 0 4)(0 0 4)(0 0 5)(0 1 5)(0 1 4)
(0 2 4)(0 3 4)(1 3 4)(1 2 4)(1 2 5)(1 3 5)(0 3 5)(0 2 5)
(0 2 6)(0 3 6)(1 3 6)(1 2 6)(1 2 7)(1 3 7)(0 3 7)(0 2 7)
(0 1 7)(0 1 6)(1 1 6)(1 1 7)(1 0 7)(1 0 6)(0 0 6)(0 0 7)

が得られ, 始めの図が描けることになった. 折角描いてはみたが, 期待したほどには美しくなかったのが残念だ.

私がHilbert曲線の話を書いたのは2009年6月19日のブログであった. そこの図で分かるように, 2次元の曲線では, 1次の曲線が縦向きに出発するなら, 2次は横向き, 3次はまた縦向き,...のように交互に出発する.

今回, 3次元のHilbert曲線を3次まで描いて分かったのは, 前回のブログの最初の図で, 赤線の1次の曲線が0からx軸方向へ出発し, 2次が0からy軸方向へ出発したのに対し, 今回の図では, 赤点から予想通りz方向へ出発していることだ. そういうものだったのだ.

上の64行にわたる512の座標は, 8行ずつ8つのブロックにして書いてある. ブロックを上から順に0,1,2,...,7ということにすると, ブロック0の座標はxもyもzも0から3である.つまり原点に近い1/8の空間を占めている. ブロック1はxが4から7なので, x方向だけが4ずれた空間にある. 最後のブロック7はzが4から7で, z方向へずれた空間にあるわけだ.

ブロック0をもう一度書いてみる.

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

赤く書いてあるのは, すべての座標が2で割り切れるもので, 各行の1つずる現れる. 各行が1つのGrayコードになっていて, その原点に近いものがそれである. 下の図の8つの赤点に対応する.

それを2で割ると

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

となり, 2次の64個の座標の1行目と同じになるが, 当たり前だ.

再帰曲線

今回は3次元のHilbert曲線が話題である. この2月は思わぬ風邪引きで, それも歳をとったせいか, 全治に時間がかかった.

それで家に引き籠もっている折, 例の「ハッカーのたのしみ」を眺めていたら, 3次元Hilbert曲線の基本部品という図を見つけた. 確かに3次元でも出来そうであるが, さてどうやったら描けるのか. しばらくそれが気になっていた.

2次元だと, 右折, 左折などを繰り返すと何となくアルゴリズムが見つかるが, 3次元の場合, これをどういえばよいか. 試行錯誤の結果, 基本部品に名前をつけ, それを組合せると描けそうだというところまで来た. もっともまだ最後まで出来たわけではなく, 今後挫折するかもしれないのだが.

とりあえず, 1次と2次の図を描いてみたのが下だ.



赤の太い線が1次, 黒の細いのが2次である. 通過した順が分かるように, 節点に0から番号を付けてある. 1次のは向こうの下の方0から始まる. この座標を000としよう. それから図で左下へ進む. この方向をxとする. つまり次の点1の座標はzyxの順に書けば001である. 今度は右下へ進む. この方向をyとする. 次の点2の座標は011だ. この伝でいくと, 010, 110, 111, 101, 100のように進むことになる.

もう1度揃えて書くと

zyx
000
001
011
010
110
111
101
100

これは誰が見てもGrayコードである.

ところで, 2次を見ると, 1次の0の回りを0→1→...→7と同様にGrayコードで進む8点で囲む. 1の回りを8→9→...→15の8点で囲む. ... 7の回りを56→57→...→63の8点で囲む. この8点の形はよく見ると5種類ある. つまり1と2(B), 3と4(C), 5と6(D)の回りの8点はそれぞれ同じで, その他に0(A)と7(E)で5種類である. Grayコードの部品がたくさんあるが, 全体がGrayコードになっているのではない.

A    B    C    D    E
zyx  zyx  zyx  zyx  zyx
000  000  011  110  101
010  001  001  111  111
110  101  000  011  011
100  100  010  010  001
101  110  110  000  000
111  111  100  001  010
011  011  101  101  110
001  010  111  100  100
x+z^ y+z^ z+xv y-zv x-zv

さてこれに分かりやすい名前を付けたい. Aを見ると最上段と最下段ではxが0から1と増えている. zは0から1になり0に戻る. yはふらふらして元に戻るから無視すると, xは増え(x+)zは010(^)の凸型ということで, x+z^と呼ぶことにする. 他のBからEも同じ要領で命名してある. 上にあった1次の形は, この方式でいうとz+y^である.

そのそれぞれを描くと次のようになる. 左上は1次の形である.



この1次の形z+y^は2次で描くとA, B, B, C, C, D, D, Eで出来ていたのだ. つまり

z+y^ = x+z^ 2(y+z^) 2(z+xv) 2(y-zv) x-zv.
        A      B       C       D     E

x座標はAで増え, Cで2回減増し, Eで減るから変化なし.
y座標はBで増え, Dで減るからy^である.
z座標はAで増減し, Cで増え, Eで減増するからz+である. 全体でz+y^となるわけだ.

2次のA, B, C, D, Eがこのように分解出来るなら, 3次も描けることになる.



分解を上の図を見ながらもう一度考えてみよう.

最初の+か-を伴う座標を主軸といおう. 次の^かvの座標を副軸といおう. 最後に名称に出てこないのを従軸といおう.

• 3と4は主軸の+, -により+か-になる. +の場合, 0,1,2では主軸は^, 5,6,7ではvになり, -なら逆になる.

• 副軸は^なら1と2で+, 5と6で-, vなら逆になる.

• 従軸の初期値は, 主軸と副軸の初期値と偶数パリティになるようにとる. それが0なら0で+, 3と4でv, 7で+. 1なら逆になる.

さっそく挑戦してみる.

x+z^ = y+x^ 2(z+x^) 2(x+yv) 2(z-xv) y-xv.

そう分かるとプログラムが書ける.

(define (xyz a) (list-ref '(x y z) a))
(define (+- b) (list-ref '(+ -) b))
(define (ud b) (list-ref '(^ v) b))

(define (bar a0 b0 a1 b1)
 (let ((a2 (- 3 a0 a1)) (b2 (modulo (+ b0 b1) 2)))
  (list
   (list (xyz a0) (+- b0) (xyz a1) (ud b1))
   (list (xyz a2) (+- b2) (xyz a0) (ud b0));0
   (list (xyz a1) (+- b1) (xyz a0) (ud b0));1,2
   (list (xyz a0) (+- b0) (xyz a2) (ud (- 1 b2)));3,4
   (list (xyz a1) (+- (- 1 b1)) (xyz a0) (ud (- 1 b0)));5,6
   (list (xyz a2) (+- (- 1 b2)) (xyz a0) (ud (- 1 b0))))));7

(z + y ^)は(bar 2 0 1 0)と入力する. すると

(bar 2 0 1 0) =>
((z + y ^) (x + z ^) (y + z ^) (z + x v) (y - z v) (x - z v))

と得られる. リストの先頭は入力パターンである. 従って上の図と同じになる.

Aを分解すると,

(bar 0 0 2 0) =>
((x + z ^) (y + x ^) (z + x ^) (x + y v) (z - x v) (y - x v))

一方, Grayコードを得るには次のプログラムでよい.

(define (foo a0 b0 a1 b1)
 (let ((b '()) (a (list 0 0 0)) (a2 (- 3 a0 a1))
  (b2 (modulo (+ b0 b1) 2)))
  (list-set! a a0 b0) (list-set! a a1 b1) (list-set! a a2 b2)
  (set! b (cons (reverse a) b))
  (list-set! a a2 (- 1 (list-ref a a2)))
  (set! b (cons (reverse a) b))
  (list-set! a a1 (- 1 (list-ref a a1)))
  (set! b (cons (reverse a) b))
  (list-set! a a2 (- 1 (list-ref a a2)))
  (set! b (cons (reverse a) b))
  (list-set! a a0 (- 1 (list-ref a a0)))
  (set! b (cons (reverse a) b))
  (list-set! a a2 (- 1 (list-ref a a2)))
  (set! b (cons (reverse a) b))
  (list-set! a a1 (- 1 (list-ref a a1)))
  (set! b (cons (reverse a) b))
  (list-set! a a2 (- 1 (list-ref a a2)))
  (set! b (cons (reverse a) b))
  (reverse b)))

初期値を設定してから, 従軸, 副軸, 従軸, 主軸, 従軸, 副軸, 従軸の順に1と0を交換する. "z+y^"は(foo 2 0 1 0)と入力する.

(foo 2 0 1 0) =>
((0 0 0)(0 0 1)(0 1 1)(0 1 0)(1 1 0)(1 1 1)(1 0 1)(1 0 0))

と得られる.

Aを１次とみて, その2次の絵を描いてみると,

なんだかうまく行っているらしい. 春から縁起がいいわ.

再帰曲線

このブログにdragon曲線のことを書いたのは, 2年以上も前であった.

dragon曲線は, 紙を同じ向きに半分半分と折り, 折り目を90度になるように広げたものであった.

最近, この折る向きを毎回逆にしたらどうなるかと思った. 紙を折ってやってみるのは, やはり3,4回が限度であり, 手元の計算機の威力を借りたくなった.

下の図を見て欲しい.



左端のAは折る前である. 下の黒丸は出発点を示す. その上の方に右向きの矢印が示すように, 最初は右に折る. するとBのようになる. 紙の長さが倍にになっているが, 今は折り方が問題なので, 分かりやすく描いてある.

出発点から途中右折したので, その角には+が付けてある.

Bを今度は左へ折る. 最初の折り目の+は左の下へ移動する. Cの出発点から辿ると, 最初は左折, 次がBで出来た右折, さらに右折する.

Cを次は右へ折るとDになる. 上の方の+や-は, 今回出来たのもで, 下の方のは左がCの下にあった折り目, 右がCの上にあったものだ.

出発点からの右折左折を抜き出すと

B +
C - + +
D + - - + + + -
この伝でつづけると次は
E - + + - - - + + - + + + - - +
となる.

一段上の列の要素を間に+と-を交互に入れたものになっているが, dragon曲線の場合と同じで, このパターンに着目する.

EをX₃, DをX₃とすると, X₃の左半分はX₂の+-を反対にしたものだ. 右半分は中央だけが違って右半分を殆んど同じである.

そこで, 'で+-の反転を表わすとすると,
X₃=X₂'+Y₂'.

結局
X₀=+, Y₀=-,
X_n=X_n-1'+Y_n-1',
Y_n=X_n-1'-Y_n-1'
となることが分かる.

Schemeでプログラムしてみる. この引数のnは上の漸化式の2^n-1になっている.

(define (r s) (if (eq? s '+) '- '+))

(define (x n i) (cond ((= i n) '+)
                      ((< i n) (r (x (/ n 2) i)))
                      ((> i n) (r (y (/ n 2) (- i n))))))
(define (y n i) (cond ((= i n) '-)
                      ((< i n) (r (x (/ n 2) i)))
                      ((> i n) (r (y (/ n 2) (- i n))))))

(map (lambda (i) (x 8 i)) (a2b 1 16))
=> (- + + - - - + + - + + + - - +)

これで右折左折の情報が得られたから, いよいよ交互折りの曲線を描くことする.

そして出来たのが次の図である. X₁からX₈までが, 色を変えながら描いてある. 左中ほど上の, 直角に下の曲がった青の線がB(X₁)である.



左から右向きに出発するのが, 出発点である. その下の緑の線がC(X₂)に対応する.

折角描いてみたが, 竜にも鳳凰にもならず, 単に三角形が出来ただけであった. 思うにdragon曲線を考えた人も, これもやっては見たが, つまらない結果だったので, このことは書いて置かなかったのかもしれない.

つまらない絵しか描けないことが分かっただけでも, 一応の知見が得られたというべきか.

再帰曲線

私が昭和13年4月に出会った算術の教科書の表紙の再帰曲線の続きである.

収束点の座標が1,1と分かったからには, さらに直観的でエレガントな求め方があるに違いないと思う.

次のように考えて見た.




A点から上向きに出発する. 4歩進みB点に至る. ここは次に下向きに出発する点である. ABの直線を赤のように引く. Bから4歩進むと, 先ほどとは上下左右が入れ替わり, C点に至るが, これは赤線の上に来るはずである. 従って, 4歩毎の点は, 右上に行ったり, 左下へ行ったりするが, 赤線の上にある. つまり, 収束するのは赤線の上である.

上向きに進み終ったD点から, 同じように考える. 4歩進むと下向きに進み終ったE点に至る. DEに青で直線を引く. Eから4歩進んだF点も同じように青線上にあるはずだ. 従って, 収束点は青線上にある.

赤線はy=x, 青線はy=1だから, 交点は1,1だ.

Gorge Polyaの「いかにして問題を解くか」には, 問題を解いた後で, もう一度振り返れと書いてある. もっとうまい解法が見つかるかも知れないからである.

とりあえずは数値計算で遮二無二解いたとしても, その後, 今回のような解法が見つかると, やはり嬉しい.

再帰曲線

はるか昔のことを思い出した. 小学校に入学したら, 算術という科目があった. (その後しばらくして, 算数と名前が変った.) その教科書の表紙に面白い図形があったのを鮮明に覚えている. 少しずつ小さくなった直角二等辺三角形が繋がっているのが見える.



(教科書の絵はウェブページwww-cc.gakushuin.ac.jp/~851051/maed/09sakurai.pdfにもあった.)

これをいつまでも続けていくとどうなるかというのが, その図形を見たときのとっさの疑問であった.

70余年を経て, 計算してみた.



これは上に伸びる三角形の列の外側の線を描いたものだ. 最初の二等辺の長さを1とする. 真上に1だけ伸びると, 45°右に折れ, √(1/2)だけ伸びるという操作を繰り返す. 8回すすむと, 次はまた真上に進むから, ここで1サイクル終了したとして, この座標を計算すると,

x座標は 0+0.5+0.5+0.25+0-0.125-0.125-0.0625=0.9375
y座標は 1+0.5+0-0.25-0.25-0.125+0-0.0625=0.9375

である. 最初に1伸びたのが, 1サイクル終ると, 次は1/16だけ真上に伸びるから, 次のサイクルの終りの座標は, xもyも 0.9375+0.9375/16 だ. さらに次はこれに0.9375/16²を足すことになるから, 初項が0.9375で公比が1/16の無限級数を計算することになる.

0.9375/(1-1/16)

だけれど, 0.9375が1-1/16だったから, この値は1だ. (上の図の赤線が1,1を示す.)

なーんだ, という結論だが, やっとせいせいした.

再帰曲線

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

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

これを X₂ + Y₂ と書くと,
X₂ = + + - + + - -
Y₂ = + + - - + - -
である. するとこの両枝は同じなので,
X₂ = X₁ + Y₁
Y₂ = X₁ - Y₁
X₁ = + + -
Y₁ = + - -
最後は
X₀ = +
Y₀ = -
この+と-のパターンはなんだろうと思っていたが, はたと思いついたのは, Gray codeであった.

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




なぜかというと, 枠で囲った部分がX₂とY₂で, その間が, 一番左の桁が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)

再帰曲線

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

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


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

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

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

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



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

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


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






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

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

再帰曲線

第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              %起動

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

再帰曲線

6月23日のblogの続き, 演習問題の解答である.

下の図の a は1枚の紙を上から見たところで, 一度も折らない, つまり0次のdragon曲線を表わす. dragon曲線には始点と終点を考える. 始点から終点へ向う矢印で, dragon曲線を表わすとすると, a は b のようになる. 矢印の先端の0は, 0次を意味する.



dragon曲線は紙を半分に折り, 折り目を直角に開くのだが, それの相当するのは, dragon曲線を表わす矢印の先に, 同じ矢印を右から合わせることである. 従って1次のdragon曲線の矢印₁は, 0次の矢印₀と矢印_0' から c のように作れる. 対応する1次のdragon曲線はd.

2次のdragon曲線と矢印の図をeに示す. 次はfのようだ.

このようにして, dragon曲線を次々と描くと, gの図が出来る. sは共通の始点である.


次にこれを青竜と赤竜で描いてみる.

aは0次の竜をつなげたもの. 青の始点は下, 終点は上で, 赤の始点は上, 終点は下である.

bは1次の青竜の終点に1次の赤竜の始点, 青竜の始点に赤竜の終点をつなげた.



0次の青を終点から1次の青の終点への上の矢印と1次の赤の終点から0次の赤の終点への下の矢印も示す. この矢印は, 前のblogの用語でいえば, 新興勢力がどこに出来たを示す.

このようにして, c, d, eはそれぞれ, 1次から2次, 2次から3次, 3次から4次へ新興勢力が拡大していく様子を示す.

これで分かるのは, 新版図は, 方向は45度ずつ時計回りにまわり, 距離は√2倍ずつ増えていることである.

一方, 例のフラクタル図の方は, 新版図は135度, 反時計回りにまわり, 距離は√2倍になっていたので, 結局同じように増殖していたことが判明した.

詳しくは, 竜の図で, 拡大したときにぶつからないことなど, 確認する必要があるが, 大体は良さそうに思った.

再帰曲線

dragon曲線と, twindragonのフラクタルの絵の関係の続きである. 「驚きだ」の理由もなんとか知りたくて, プログラムをいろいろ修正し, ヒントになりそうな絵を描いてみた. 例えばこれ.



フラクタルの図には, 黒丸の他に白丸がある. dragon曲線では, 太線と2重線がある. 想像出来るように, 白丸は0, 1, 2,...と点をプロットした順で, 後半のものである. また2重線も, dragon曲線を順に描いた時の, 後半の線である. つまり白っぽいものは, いわば新興勢力である.

フラクタルの新興勢力は, i-1の冪乗の地域に発展する. それに対し, dragon曲線でも, 顕著な向きは分からぬまでも, 大雑把な方向は読める. すると左上の次数n=1から順に

n  赤の新勢力 青の新勢力 フラクタルの新勢力
1       右          左            右
2      右下        左上          左上
3       下          上            下
4      左下        右上          右上
5       左          右            左

となり, フラクタルの新興勢力は, 赤か青の新興勢力の方向と合っている. 当然青と赤の新興勢力は, 逆向きだ.

dragon曲線の新興勢力が, どういう向きに出来るかは未検討だが, 「低気圧が来たので天気が崩れる」程度の, 一応の説明にはなったかもしれない.

しかしこれは, 「高気圧は下にあり, 台風は右上に進む」といっているようなもので, もう少し定量的な説明は出来ないかと思った. そして注目したのは, 折り目の点である. dragon曲線は, もともと紙を半分に, 半分に,...と折って作ったことになっているので, その辺に再帰のヒントはないかと考えた.

下の図は64辺からなるdragon曲線で, 端点と曲がり角には0から順に番号を振ってある. 最初の折り目は, 左中央やや下の32番である. その次は16番と48番になる. そこで32番から16番と48番への関係をみると, ちょうど真上と真横にあった(赤線). 次に16番から8番と24番を見ると, ちょうど斜め方向にあった(緑線). 次を見ると再び真下と真横(橙線), さらに斜め線が見えてきた(青色). しかも隣りとの長さは, 順に4, 2√2, 2, √2である. これでもうtwindragonの仲間であることがほとんど判明した.



図がごちゃごちゃしないように省いた, 最後の折り目間の関係は, 点4k+2からそのプラスマイナス1へである. つまり2から1へと3へ, 6から5へと7へと,...,62から61へと63へと, である. 描いてはないが想像出来る.

ここまでくると, TAOCP風の筆法に従えば, 後は「演習問題参照」となる.

それにしても, こういう絵がすぐ描けて, テストも出来るし, 説明も出来るし, お絵書きプログラムは, 私にとって神様仏様だ.

再帰曲線

dragon曲線と, twindragonのフラクタルの絵は, 関係がありそうでなさそうである. TAOCPには "B. Mandelbrot named S the "twindragon" because he noticed that it is essentially obtained by joining two "dragon curves" belly-to-belly." と書いてあるが(演習問題4.1-18の解答), その意味は, 前回引用したWilliam Gilbertの論文を見て分かった.

前回のtwindragonの絵は, i-1進法の小数であったが, i-1進法の整数でも話はほとんど同じである.(...a₃a₂a₁a₀)₂は

n = &Sigma_k=0 a_k*(i-1)^k

を表わす.


整数版のtwindragonの絵を描くには, まず原点に...=a₃=a₂=a₁=a₀=0に対する0の点を置く. 次にa₀だけが1で他が0の1の点は1*(i-1)⁰=1 なので, (1,0)に置く.

前回にも書いたように, 次の2と3は, 0と1の点を(-1,1)だけ移動するのだが, これは(i-1)¹=i-1だからである. 上の図の1の矢印の向きである.

冪が2になると, これは1の矢印の2乗で, (i-1)²=i²+2i+1=-2iの方向へ, 0,1,2,3の点を移動する.

i-1の冪乗は, いちいち計算することはない. 昔々習った複素数の乗算z₀*z₁=zでは, zの偏角は, z₀の偏角とz₁の偏角の和だし, zの絶対値は, z₀の絶対値とz₁の絶対値の積であった.

だから, 2の矢印は, 1の矢印をもう135度回して下向きにし, 1の矢印の絶対値√2の2乗で2の長さにする. それが矢印2である.

矢印2に矢印1を掛けると, 偏角はまた同じだけまわり, 右上45度を向く. 絶対値は2の√2倍だから, 2掛ける2の正方形の対角線になる. 矢印3の先は, 点(2,2)である.

要するに次々の矢印は, 135度回転し, 長さを√2倍すればよいことが分かる.

√2倍というのは, A版, B版の紙の大きさを連想させる. つまりAn版の紙を2枚横に並べる. その上にAn-1版の紙を並べる. その横にAn-2版の紙を並べる. ... twindragonの絵はこういうものだったかと考えると, 有難みは減る. そういえば規格の紙の寸法は, フラクタルの身近な例であったのだ.



さて, dragon曲線とフラクタル図の関係である. 上の図に, モンゴル出身力士のような, 青竜と赤竜を示す. 点対称になっている. またフラクタル図も示した. 青竜と赤竜の端点はちょうど左右に並ぶように描いてある. これを端点が繋がるようにずらし, さらにフラクタルもずらすと, 竜の縦線とフラクタルの点が一致するのである. 驚きだ!

ついでにいくつかの竜についても描いたのが, 下の図である.



こういうことに気づくとは, さすがMandelbrotである.


1958年秋, IBMのOssiningの研究所に, 先輩の蒲生秀也先生を訪問した時, 先生はMandelbrotがその研究所にいるので, 紹介しようといわれたが, その日Mandelbrotはたまたま休暇で, 会うことは出来なかった.

再帰曲線

dragon曲線ですぐ思い出すのは, TAOCP(The Art of Computer Programming)4.1にあるtwindragonのフラクタルの絵である.



右上に灰色の竜がいる. その左下に黒の半分の大きさの竜がいる. さらにその下に灰色の半分の大きさの竜がいる. この後, 黒, 灰, と交互に小さくなり, 最後は同じ大きさの黒と灰が並んで終る.

よくよく見ると, 一番大きい灰色の竜は, 相似形の小さい竜で出来ている. 小さい竜は128個ある. 黒竜はよく見えないがこれも同様な小竜で出来ている. 最後の竜は, 小竜1個だ.

種明かしをすると, この絵は i-1進法の小数を示すものだ. 通常の2進小数(.a₁a₂a₃a₄...) ₂(a_i=0,1)は,

n = &Sigma _i=1 a_i*2^-i

を表わすが, i-1進法では, この式の2の代りにi-1を使う. iはいうまでもなく虚数単位である.

これらの数の表わす値は,

0.1 =-1/2-1/2i
0.01=+1/2i
0.001=1/4-1/4i
0.0001=-1/4
0.00001=1/8+1/8i
0.000001=-1/8i
0.0000001=-1/16+1/16i
0.00000001=1/16

で, 従って

0    =0
0.00000001=1/16
0.00000010=-1/16+1/16i
0.00000011=+1/16i
0.00000100=-1/8i
0.00000101=1/16-1/8i
0.00000110=-1/16-1/16i
0.00000111=-1/16i

これを複素平面に描くと下のようになる. これは面白いことに, 複素平面を重複せずに覆う.



構造的にはまず中央に0の点が出来る. 0.00000001=1/16なので, この0の点を右へ1単位移動すると1になる. 次に 0.0000001=-1/16+1/16iなので, 0と1の点を左へ1, 上へ1移動すると2と3になる. また0.000001=-1/8iなので, 0,1,2,3を下へ2単位移動して, 4,5,6,7とする. このように, それまで出来た点を次々と移動すればよい. 0.1=-1/2-1/2iなので, 黒い竜は灰色の竜の左下にいたわけだ.

これを256点までとると, 下のようになる.



最初のフラクタルの絵で, 灰色と黒の小竜が並んでいたのは, この図では下から6段目, 一番虚数軸に近い2個である.

Schemeは, 複素数は自家薬籠中の機能なので, こういう計算はなんでもない.

(define mydevice (make-graphics-device 'x))
(define (n2b n m)
 (if (= m 0) '()
     (cons (modulo n 2) (n2b (quotient n 2) (- m 1)))))
(define (evaluate n)
 (if (null? n) 0
   (+ (/ (car n) b) (/ (evaluate (cdr n)) b))))
(define b -1+i)

(do ((i 0 (+ i 1))) ((= i 256))
 (let* ((a (reverse (n2b i 8))) (p (evaluate a))
        (x (real-part p)) (y (imag-part p)))
   (graphics-operation mydevice 'fill-circle x y 0.02)))

(graphics-draw-line mydevice -1 0 1 0)
(graphics-draw-line mydevice 0 -1 0 1)

もうこれはM.C.Escherの世界でもある. 詳しい話はwww.math.uwaterloo.ca/~wgilbert/Research/MathIntel.pdfをご覧あれ.

再帰曲線

普段使うプログラミング言語はSchemeが多い. その1つの実装, MITSchemeにはgraphicsの機能があることは, うすうす知ってはいたが, 使ってみたら存外簡単であった.

(define mydevice (make-graphics-device 'x))

で四角い描画領域が出来る. 座標は左下が(-1,-1), 右上が(1,1)だ.

点(x0,y0)から点(x1,y1)へ直線を引くには,

(graphics-draw-line mydevice x0 y0 x1 y1)

でよい.

早速絵を描いてみる. 関数型言語だから, 再帰的な絵がいい. となるとhilbert曲線とdragon曲線である. 描き方は後回しとして, 出来映えは





の通り.

どちらも良く知られている曲線であるが, 私は次のように考えて描いた. hilbertの方は, 以前, 情報処理学会誌にも書いたことがある.



上の図で, 左端は, 1次から4次までのhilbert曲線が重ねて描いてある. こうして見ると, フラクタルのように, 最初の直線を, 少しずつ曲げて作られていることが分かる. 中のは, hilbert曲線の要素で, 見かけは同じようだが, 下のXは左から来て, 左折, 右折, 右折, 左折し右に抜け, 上のYは右から来て, 右折, 左折, 左折, 右折し左へ抜ける.

中の下のXを右のXのようにするには, +を左折, -を右折とすると,
X -> + Y F - X F X - F Y +
Y -> - X F _ Y F Y - F X -
と変換すればよい. Fは1歩前へ進むことである.

従って, hilbert曲線のプログラムは, 関数名を変更し, +をleft, -をright, Xをp, Yをqと書くと,

(define mydevice (make-graphics-device 'x))

(define (l) (let ((t dx)) (set! dx (- dy)) (set! dy t)))
(define (r) (let ((t dy)) (set! dy (- dx)) (set! dx t)))
(define (f)
 (graphics-draw-line mydevice x y (+ x dx) (+ y dy))
 (set! x (+ x dx)) (set! y (+ y dy)))

(define (p) (if (> n 0) (begin (set! n (- n 1)) (l) (q) (f)
 (r) (p) (f) (p) (r) (f) (q) (l) (set! n (+ n 1)))))
(define (q) (if (> n 0) (begin (set! n (- n 1)) (r) (p) (f)
 (l) (q) (f) (q) (l) (f) (p) (r) (set! n (+ n 1)))))

(define n 6) (define x -0.8) (define y -0.8)
(define dx 0.025) (define dy 0)
(p)

となる.

次にdragon曲線.



dragon曲線は, 紙を半分に折り, それを開くと, 上の左端のようになる. 次に紙を半分に折り, さらに半分に折り, それを開くと, 左から2番目のようになる.


線の開始点には, 小さい黒四角, 折り返し点には小さい白四角がついている.

半分に折る方向は, 山折, 谷折両方があるが, ここでは開いた絵が, 左折するようにしている.
そして左折を+, 右折を-と表示する. 折り返し点は常に + にしてある.

半分折りを3回続け, それを開くと, 右端のようになる. この辺からがdragon曲線らしくなる. 曲がる向きを順に書くと

+ + - + + - - + + + - - + - -

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

つまりこれを X₂ + Y₂ と書くと,
X₂ = + + - + + - -
Y₂ = + + - - + - -
である. するとこの両枝は同じなので,
X₂ = X₁ + Y₁
Y₂ = X₁ - Y₁
X₁ = + + -
Y₁ = + - -
最後は
X₀ = +
Y₀ = -

これをプログラムにすればよいから, 前と同様に, +をleft, -をright, Xをp, Yをqと書くが,

(define mydevice (make-graphics-device 'x))

(define (l) (let ((t dx)) (set! dx (- dy)) (set! dy t)))
(define (r) (let ((t dy)) (set! dy (- dx)) (set! dx t)))
(define (f)
 (graphics-draw-line mydevice x y (+ x dx) (+ y dy))
 (set! x (+ x dx)) (set! y (+ y dy)))

(define (p) (if (= n 0) (begin (l) (f))
 (begin (set! n (- n 1)) (p) (l) (f) (q) (set! n (+ n 1)))))
(define (q) (if (= n 0) (begin (r) (f))
 (begin (set! n (- n 1)) (p) (r) (f) (q) (set! n (+ n 1)))))

(define n 10) (define x 0.4) (define y -0.2) 
(define dx 0.025) (define dy 0)
(f) (p)

となる.