2009年6月23日火曜日

再帰曲線

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風の筆法に従えば, 後は「演習問題参照」となる.

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

2009年6月22日月曜日

再帰曲線

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進法の整数でも話はほとんど同じである.(...a3a2a1a0)2

n = &Sigmak=0 ak*(i-1)k

を表わす.


整数版のtwindragonの絵を描くには, まず原点に...=a3=a2=a1=a0=0に対する0の点を置く. 次にa0だけが1で他が0の1の点は1*(i-1)0=1 なので, (1,0)に置く.

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

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

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

だから, 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はたまたま休暇で, 会うことは出来なかった.

2009年6月20日土曜日

再帰曲線

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



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

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

種明かしをすると, この絵は i-1進法の小数を示すものだ. 通常の2進小数(.a1a2a3a4...) 2(ai=0,1)は,

n = &Sigma i=1 ai*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をご覧あれ.

2009年6月19日金曜日

再帰曲線

普段使うプログラミング言語は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個あるが, 中央の + が最初の折り目で, その両側はもともと重なっていたのを開いたから, 並び順が中央を中心にして対象で, しかも+と-が入れ替わっている.

つまりこれを X2 + Y2 と書くと,
X2 = + + - + + - -
Y2 = + + - - + - -
である. するとこの両枝は同じなので,
X2 = X1 + Y1
Y2 = X1 - Y1
X1 = + + -
Y1 = + - -
最後は
X0 = +
Y0 = -

これをプログラムにすればよいから, 前と同様に, +を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)

となる.

2009年6月15日月曜日

手回し機械式計算機

手回し機械式計算機の話である. 手回し機械式計算機を知らない人も, 使ったことがない人も多いと思うので, まずその簡単な説明から.



この写真は東京理科大学の近代科学資料館で撮らせていただいた, 非常に古いBrunsvigaである. 上の筒状の右が置数レジスタ(S)で, レバーを動かして数値を置く. 手前の横長の部分(キャリッジという)には, 左に回転数レジスタ(M), 右に結果レジスタ(R)がある. 写真では見えないが, 筒状の向うにクランクハンドルがあり, 正方向に回転すると, 置数レジスタの数が結果レジスタに足され, 回転数レジスタが1増える. 負方向の回転では, 置数レジスタの数が引かれ, 回転数が1減る. 乗算をするには, 置数レジスタの数を, 1の桁, 10の桁, ...に足すので, キャリッジは左右に移動できる.

123に456を掛ける場合, 置数レジスタに456を置き, 結果レジスタを0にし, 1の桁に3回, 10の桁に2回, 100の桁に1回足す. 結果56088が得られる. これに更に789を掛ける場合, この古い機械では出来ないが, 後年の機械では結果レジスタの値を置数レジスタに移すback transferという機能があり, 56088が置数レジスタに置ける. 1000の桁に1回足し(1000倍が出来た), 100の桁で2回引き(800倍), 10の桁で1回引き(790倍), 1の桁で1回引く(789倍). これをshort cut乗算という

昔はこういうことをやっていたのもだ. 開平も出来ることは, 2009年1月21日のブログに書いた.

ところで, Brunsvigaには, もっと多機能な機械があった. それの使い方を考えたい.

その計算機はBrunsviga DuplaBrunsviga Doppelである.

これらの計算に1930年前後に登場した. 最初の写真の計算機と違うのは, 回転数レジスタが置数レジスタの上に来た; 置数レジスタの数値が良く見えるように, チェックレジスタが置数レバーの上に新設されたことである.

さてDuplaは, 置数レジスタと結果レジスタの間に, 第2結果レジスタ(R2)を置いた. ハンドルを回転すると, 置数レジスタの数値は, 第1結果レジスタ(R1)には加減されるが, 第2結果レジスタには, 加減算を止めることが出来る. また, どちらの結果レジスタからも, back transferが出来る. さらに第1結果レジスタの下に, thumb wheelという仕掛けがあり, 第1結果レジスタに, 直接数値を入れることが出来る. この機構は除算の時に便利で, これがないと, 非除数をR1に置くのに, 置数レジスタに一旦被除数を置き, 1回加算しなければならず, とくに長い被除数だとその加算が2回になる.

Duplaの使い方については, Comrieの書いた論文がウェブで見つかる. それによると, 第2階差まで使った数表が作れるというので, やってみよう.

y=x2+x+41の表と, その第1階差, 第2階差を示す. (これは素数生成(!!)の2次式である.)


x y Δ' Δ''
0 41
2
1 43 2
4
2 47 2
6
3 53 2
8
4 61 2
10
5 71 2
12
6 83 2
14
7 97 2
16
8 113 2
18
9 131

Comrieの方法によると,
初期設定
関数値41をR2に置く. 第1階差2をSに置く. 第2階差2をR1に置く. Mを0にする.

ステップの進め方
ハンドルを1回転する. するとR2に新しい関数が出来る. Mに引数が出来る. R1に第1階差が出来る.
従って, そこでR2を記録する. R1をSにback transferする. R1に第2階差2を置く.


次々と進めると

x y Δ' Δ''
1 43 4 2
2 47 6 2
. .. . .


となり, 上の表が出来る. R1に定数の第2階差を毎回置くのは面倒なような気もするが, 第2階差が規則的に変動しても計算に使えるので, これは意味がある.

一方, Doppelは大雑把にいうと, 計算機を2台, 横に連結したものである. クランクハンドルは1個だが, 回転数レジスタ, 置数レジスタ, 結果レジスタは左右に1個ずつあり, クランクハンドルを正方向に回転すると, 普通には, 右置数レジスタの値は右結果レジスタに足され, 左置数レジスタの値は左結果レジスタに足される. 負方向に回転すれば, 足す代りに引かれる. 中央に上向き矢印が2個並んでいる絵と, 上向きと下向きが並んでいる絵があり, そのクラッチを切り替えると右では足し, 左では引ける. さらに接続を切ることも出来る. またDuplaのように, 結果レジスタに直接設定も出来る.

では, こういう重連計算機はどういう時に使うか. あまり使い方を書いたものは知らないが, 1つの利用法を考えた.

十進法の1023を八進法ではどうなるか. もちろん1777である. 以下がBrunsviga Doppelによる計算法である. 左と右の回転は逆にする, つまりハンドルは正方向に回し, 右は加算, 左は減算する.


行 左S 左R 右R 右S
0 512 1023 0 1000
1 64 511 1000 100
2 447 1100
. .. ... .... ...
7 127 1600
8 8 63 1700 10
9 55 1710
. . .. .... ..
14 15 1760
15 1 7 1770 1
16 6 1771
.. . . .... .
21 1 1776
22 0 1777


とまぁこういう使い方もあるわけだ. もう少し詳しくいうと, 左置数レジスタ83, 右置数レジスタに103を置く. 左Rに1023を置き, 右Rを0にする. 左Rから左Sが引ける限りSを引く. 引けなくなったら(1行目), 左Rを82, 右Rを102にして, また引けなくなる(8行目)まで引く. これを左Rが0になるまで繰り返すと, 右Rに八進表現が出来ている.