2009年12月29日火曜日

微分解析機

微分解析機(differential analyzer)は変数の値を軸の廻転で表す. 変数同士の和は, 差動歯車(differential gear)を使うといわれている. 英語はおなじdifferentialでも意味は違う.

差動歯車は自動車の駆動装置に使われていることでよく知られている. 息子が小学校の上級のころ, プラモデルで無線操縦のレーシングカーを作るのを見ていたら, プラスティック製の差動歯車の部品があり, それで差動歯車を組み立て, 中にグリースを詰めて密封していた. こんなところまでこだわるのかと恐れ入った.

差動歯車の構造を下のようである.

四角で表す箱状のものに, かさ歯車が3つ入っている. 左右のかさ歯車はそれぞれ外部のx軸, y軸に連結している. 図に見るように, x軸とy軸は直接にはつながってはいない. 中央下のかさ歯車は, 左右のかさ歯車に噛み合い, 軸は箱状のものに刺っている. x軸, y軸が同方向に同速度で廻転すると, 中央のかさ歯車は, 自分の軸の周りには廻転せず, 箱状のものをx, y軸の周りに, 同速度で廻転させる. 箱状のものがx, y軸の周りに廻転すると, それと一体の歯車もxやyと同様に廻る.

xの廻転を止め, yだけ廻転させると, 下のかさ歯車が廻転しながら, 箱状のもの(z軸といおう)も廻転する. しかし今度はyの半分である. z=(x+y)/2なのだ.

作動歯車の動きを滑車で表すと次の図のaのようになる. 滑車に紐をかけ, 両端をxとyだけ引くと, 滑車(下のかさ歯車に相当)は廻転せずに同じ長さだけ下がる. bはxを引かない時で, 滑車は廻転し, yの半分だけ下がる.


滑車が廻転するのは, 紐が滑らないことなので, 紐の代りに滑らない棒を使っても同じなことをcは示す. これはxとyのかさ歯車を展開したものである.

図では箱状の歯車を半分の歯の歯車に噛ませて(2倍して), x+yが出来ることを示す.

自動車のいわゆるデフは, このx+yの軸が駆動軸で, 両車輪xとyを駆動する. エンジン(ミッション)からの駆動軸は車軸と直交するので, この歯車はかさ歯車である.

ところで通常デフといわれる装置は, 最初に書いたように, ディファレンシャルである. ディジタルをデジタルというのに似た変形である.

佐々木達治郎先生の「計算機械」(1947年刊)に東大航空研で試作した微分解析機の構造が説明してあり, その中に機素の加算に遊星歯車(planetary gear)を使うと書いてあった.

遊星歯車といえば, 内装式の自転車の変速機が有名である. 昔の通信機器は, 同調ダイアルにやはり遊星歯車を使っていた. つまり背後に丸い文字版があり, 手前に黒いつまみがある. つまみを廻しても, 文字版ははるかに遅く廻り, 微妙な同調が出来るようになっていた. 文字版には副尺もついていて, 高級感を出していた.

下が遊星歯車の図解である.

中央の太陽歯車の廻りに同じ大きさの遊星歯車があり, 遊星歯車は太陽と外側の(内向きに歯を持つ)歯車に噛み合う. また遊星歯車の軸は遊星キャリヤに乗っていて, 歯車の位置が変るとキャリヤの軸も廻転する. (こういう絵を描くのは楽しい.)

これでどうして差動歯車になるか, 即座には分からないが, もう1つの絵を描いてみると納得出来る.


これは太陽歯車と, 外側の歯車を1列に延したもので, 一応遊星歯車と等価を考えられる. 太陽歯車の廻転は下の棒状の歯車(ラック)の移動, 外側のは上の歯車の移動であり, キャリヤの軸の廻転は中央の歯車の軸の左右の移動と思えばよい.

歯車の場合は, 廻転角度が重要な情報だが, ラックでは移動の歯数で考えることになる.

こういう絵をみると, 上の差動歯車と同じになり, なるほど加算機構に使えるを分かる.

差動歯車の解説を読むと, かさ歯車のものより, 遊星歯車型の方が基本のように書いてあった. われわれは車の知識があるので, かさ歯車が一般的を思い込んでいたようである.

2009年12月21日月曜日

微分解析機

機械式の微分解析機では, 変数同士の乗算は, 積分機を2台と加算でするのだと思っていた. ところが図面に乗算卓というのがあり, 積分をせずとも乗算しているらしかった.

幾何学での乗算法はよく知られている.



点Pから半直線を引き, PA=a, PB=bとなる点A, Bをその上に取る. またPから別の半直線を引き, PC=1の点Cを取る. A,B,Cを通る円Oを描き, PCの延長上で円との交点をDとするとPDの長さxがa*bになっている.

これは方冪の定理による. その説明が次の図である.



PO2-AO2を考える. 正三角形POMでPythagorasの定理を使うと
PO2=PM2+OM2
同様に
AO2=AM2+OM2
辺々相引くと
PO2-AO2=PM2-AM2=(PM-AM)*(PM+AM)=PA*PB
つまりまえの図のPA*PB=a*b=PO2-AO2
=PO2-半径2=一定

である.

しかし微分解析機でこれをやっているか疑問であった. ところである文献に微分解析機の使い方の説明があり, 乗算卓の記号が下の図のようになっていた.



元々の記号に赤や青の色はない.

微分解析機からs, t, xの3本の軸で接続されている, s*t=xと出力が出るらしい.

x軸の右にはハンドルのような絵がある. そこで考えてみて, こうなっていると推測した. t軸が廻転すると, 軸に接続されている赤色のT型の足が円の中心を軸に振れる. 一方s軸の廻転で青い棒が左右に移動する. オペレータはハンドルを廻し, 青軸上のカーソルを赤線の上に保つのである. そのハンドルの廻転がx軸で出力される.

その説明が次の図だ.



円の中心がAの点である. tもsもAから左右に計測する. tが増えるとB点は左へ移動し, ABと直角についている腕, ADも一緒に廻転する. EDの長さをxとすると, 三角形ABCと三角形ADEは相似であり, ACの長さを1とすると, 1/t=s/x 従ってx=s*t である.

そういう乗算卓を使っていたとは面白い.

2009年12月12日土曜日

微分解析機

その昔, 微分解析機という低速のアナログ計算機があった. 1階の常微分方程式を解く機械式の積分機を何台かと, 入出力テーブルを持ち, その間の情報は多くの廻転軸で伝達した.

Wikipediaによれば, 原理を発明したのは1876年, James Thomsonで, 実用機の製作は1927年, MITのVannever Bushによる.

我が国にも2, 3台はあったらしい. 神楽坂にある東京理科大学の近代科学資料館にいくと, 展示されている. 私は1950年代に, 東大生産技術研究所にあった大型機を見学したことがある.

ないよりは数段ましの計算機だが, 大阪大学にいらした, 清水辰次郎先生(1887-1992)が微分解析機の使用経験を書かれたものを読むと, これを使うのは, 結構大変であったらしい.

清水先生の話は

d2x/dt2-μ(1-x2)dx/dt+x=0

を積分機2台の微分解析機でどう解いたかというものである. この微分方程式はVan der Polの式といって, 非線形屋さんの好きなものである.

2階の式は2つの変数に分け, 1階の連立にするのが常道である.
dx/dt=y      積分機0番
dy/dt=g(x,y)    積分機1番
where g(x,y)=μ(1-x2)y-x

xの初期値を1, yのを0, μを0.1にして, この結果をxとyとの図にしたのを下に示す.


今ではPCで簡単に図が描けるが, 昔は手間がかかったので, こういうきれいな図が得られると, 額にいれた飾っていた人もいた.


さて, 阪大の微分解析機は積分機が2台あったから, これで解けそうに見えるが, 問題はμ(1-x2)y-xの計算である. 微分解析機の変数は前述のように廻転軸なので, 最後の引き算は差分ギアを使えばよいし, 最初のμの掛け算はμが定数なので, ギア比で実行できる. xの自乗やそれにyを掛けるように, 変数同士の乗算があると, 微分解析機は積の微分の式を使い, 2台の積分機を使うことになるのだが, 阪大には余っている積分機はない.

そこで聞くも涙の話になるのだが, 清水先生は, 固定したμと, ある範囲のxとyとに対して, g(x,y)をあらかじめ計算し, 表か図にした.

適当な初期値から, (例えばx=1, y=0から,) 起動するとき, それに対応するg(x,y)を積分機1番に入力する. するとx, yが計算され, 出力テーブルにプロットされる; その値を読んで1人が叫ぶと, もう1人がそれに見合うg(x,y) を入力するという, 恐ろしいことをやった. この遅れが誤差になるのは承知の上だが, 傾向をみるだけにはそうするしかなかったのであろう.

私の興味は, g(x,y)はどういうものかにあった. 微分解析機では, f(x)のような1変数の入力なら, 入力テーブルにこの関数形を描いておき, 人手でトレースするなり, 生研の機械なら, 光電管で曲線を自動追尾するなりして, 入力出来るが, gは2変数関数なので, 3次元の入力装置を作れば自動化できたかも知れない. それで関数の形が知りたかったのである.

-4≤x,y≤4について, 0.2おきにg(x,y)を描いたのが下である. もちろん3次元には表示できないので, 高さを向こうへ倒して表現してある. つまりy=4は一番上の線で, 基準が赤線であり, x=-4のとき, g(-4,4)=-2と読む. 他の3つの隅の例では, g(-4,-4)=10, g(4,4)=-10, g(4,-4)=2である.


これならなんとか立体模型をつくり, 入力装置も出来たのではないかと想像する.

この後, この図を等高線で描いてみた. 等高線のプログラムは, 天気図を描きたいという人などが, 昔から挑戦しているが, 私は始めてプログラムしてみた.

それが次の図である.


さらに, 横倒しの図と等高線図を重ねたものも示す.



積分だけでなく, こういう描画も楽になったとつくづく思う.

2009年11月29日日曜日

彩色立方体

東急ハンズで1辺が2センチの木の立方体を売っていたので, 早速買い求めた. 50個で500円程度であった.

これに白, 赤, 緑, 橙, 青, 黒の紙を張って30個の彩色立方体をこしらえた.



作り方は次の通り.

まず30個の立方体のそれぞれの各面に, 色の名を鉛筆でつける.
白で1辺が3センチ程度のほぼ正方形の紙を30枚用意する.
糊で立方体の白の面にその紙を張る.
糊が乾いたら, カッターで余白を切り取る.
白の作業が終わったら, 他の色でも同様に行う.

と書けば簡単だが, 180枚についてこの作業をするのは, 一仕事であった.

今回の反省は, やはり緑, 青など色が濃すぎた; 糊が隣りの面にはみ出て汚れたところがあった; カッターの刃が傾き, 稜を削った個所があった; など.

プログラムの結果と合わせるための, 立方体の番号を探すカタログを用意する必要も感じている.

2009年11月21日土曜日

彩色立方体

Martin Gardnerの本に30 color cubesという話題があった.

立方体には面が6つある. それに相異なる6色を塗る. 廻転して同じになるようなものは作らないとすると, 30通りの塗り方がある.

6色を0,1,2,...,5とする. 立方体の面にも名前をつける. 島内剛一先生の「ルービック・キューブと数学パズル」にあったように, 上をtopのt, 下をbottomのb, 北をn, 東をe, 南をs, 西をwとするのは, 常識だろう. これを(t n e s w b)の順に並べることにする.

まずtを基準として色0を割当てる. するとbには, 他の5色のいずれかが来る. そして残りの4色で側面を塗るが, 4色の最小の色をnに割当てると, 後は3色の置換の6通りになり, 都合30通りというわけだ.

こんなプログラムを書いた.

(define (gencubes)
(let ((t 0) (neswb '(1 2 3 4 5)))
(apply append (map (lambda (b)
(let* ((nesw (delete b neswb)) (n (apply min nesw))
(esw (delete n nesw)))
(map (lambda (esw) (append (list t n) esw (list b)))
(permutation esw)))) neswb))))

2行目でtを0, neswbを(1 2 3 4 5)とし, 3行目のmapのlambda変数でそのいずれかをbとする. neswbからbを除いたものがneswで, その最小をn, 残りをeswとした. それの辞書式順の置換を作り, 前に(t n), 後に(b)をappendして30通りを作っている.

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

これを辞書式順にソートしたものを, 0は白, 1は赤, 2は緑, 3は青, 4は橙, 5は黒として表示したのが, 以下の図である.

左上の0番で説明すると, それぞれは2つの6角形が連接した形になっていて, 左の6角形は立方体を北東の上方から眺めたもので, t,n,eの3つの面がそれぞれ菱形で表示してある. 上にt=0, 右下にn=1, 左下にe=2の面が見える. 一方, 右の6角形は,南西の下方から眺めたもので, 右上にs=3, 左上にw=4, 下にb=5の面がある. wの4とnの1は, 図でつながっているだけでなく, 実物でもつながっている. 立方体の側面は, 1,2,3,4と北から右まわりにつながっている様子が描いてある.

そこでクイズなのだが, 仮にいま話題にした0番を見本とすると, 他の立方体から8個を選び, 2x2x2のちょうど2倍の立方体を作り, その立方体の外方の面を見本と同じ色にし, 立方体の内側で相接している面同士は, 同じ色にするというのである.

通常なら, 積み木を用意し, とっかえひっかえやってみることになるわけだが, そこはやはりプログラムで全解探索せざるべからずである.

そしてやってみた結果が下の図で, 解は2通りあった.

この図の見方は, それぞれの解で, 大型の菱形に配置した4個の立方体が上下2段に描いてある. 上の4個がそのまま上段に, 下の4個が下段になる. それぞれの菱形の下(0と4)が北東, 左(1と5)が南東, 上(2と6)が南西, 右(3と7)が北西を占める. それぞれの解で, 上段の4個(0,1,2,3)のtは0, 下段の4個(4,5,6,7)のbは5, 北の4個(0,3,7,4)のnは1,...のように, 外の6面は4個ずつ見本と同じである. また相接する面も, 例えば上の解の1のb=4は5のt=4で一致しているし, 1のn=5は0のs=5で一致している.

プログラムは通常の深さ優先の探索木を辿るものである. 解の図のように, 8個の立方体を0,1,2,...の順で候補を探す. そのため, ある立方体をあらゆる方向から見た24個を作るプログラムと, (0 1 2 () () ())のようなパターンと, ある立方体の色の順が, 合っているかを見るプログラムを用意した.

探し方は以下の通り.

探すパターン 決った色
c0 (t n e - - -) s0 w0 b0
c1 (t s0 e s - -) w1 b1
c2 (t - w1 s w -) n2 b2
c3 (t n w0 n2 w -) b3
c4 (b0 n e - - b) s4 w4
c5 (b1 s4 e s - b) w5
c6 (b2 - w5 s w b) n6
c7 (b3 n w4 n6 w b)


30通りの立体の組を全候補とする. その中から見本を選ぶ. 上の例では0番. 見本からt, n, e, s, w, bを設定する. 全候補から見本を除いたものをc0用候補とし, そこからc0のパターンに合ったものをc0とする. c0が決るとs0, w0, b0も決る. c0用候補からからc0を除いたものをc1用候補とし, そこからc1のパターンに合ったものをc1とする. c1が決るとw1, b1も決る. 以下同様にしてc7まで決ればc0,c1,...c7を出力する.

プログラムはこれだけである.

30種の立方体は, 全て対称的に出来ているから, 0番を見本とする解が得られれば, 1番,...,29番を見本とする解もまったく同様に作れる. 実際プログラムがあるのだから, 全部やってみた. 全く同じような解が2つずつ得られた.

追記 上と下の解で使われている立方体の番号(上の図参照)は, それぞれ

6 4
22 7 1 16
19 12
12 19
16 1 7 22
4 6

であり, 2つの解で同じ組み合せであった.

2009年11月8日日曜日

二進木の置換表現

図のような二進木があるとする. ノードの番号は中順序で付けてある.

二進木だから左リンクと右リンクがある. つまり

k 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
LLINK[k] 0 0 1 0 0 4 0 0 8 0 7 6 0 13 0
RLINK[k] 2 0 12 5 0 11 9 0 10 0 0 14 0 15 0

この2組のリンクを1組に圧縮する方法がTAOCP V4F4にあった. permutation representation of binary treeという.

ノードがn個あれば, 左と右でリンクは2nあるわけだが, 意味有るリンクは1を子にするもの, 2を子にするもの, (3はルートだから, 子にするものはない,)4を子にすもの,... という n-1個しかない. あとは子なしのnilのリンクだ. n=15の上の表でも, 0 は16ある.

この0は無駄なので, それを使おうと考えるのが, この方法の味噌である.

ノード番号は中順序なので, LLINK[k]<k かつ LLINK[k]!=0ならLLINK[k]は重要な情報である.

LLINK=0の場所がうまくRLINKに転用出来るかが, 問題である. つまり図で,
1, 3, 4, 6, 7, 9, 12, 14 (1)
から出ている右リンクの情報を,
1, 2, 4, 5, 7, 8, 10, 13, 15 (2)
の左リンクの場所を使って表現出来るかということだ.

ありがたいことに, 右部分木には, 必ずもっとも左のノードがあり, もっとも左だから, 左リンクは空いていて, しかも中順序だと, ノード番号が1だけ多い. つまり上の(2)の列は, 先頭の1を除いて(2)の列より1ずつ多い.

例えはノード3の右部分木は12をルートとするものだが, そのもっとも左のノードは4. ノード4の右部分木のもっとも左のノードは5. ノード6の右部分木は11をルートとするものだが, そのもっとも左のノードは7である.

これを活用することにすると, 3の右リンク12を4に置かせて貰うが, 部分木のルート自身がもっとも左でなければ, いまの場合がそうだが, 4<12になり, 4の右リンク5を5に置かせて貰う場合, 右部分木のルートがもっとも左なので, 5≤5となる.

要するに, こういう情報だけを保存すると,

k 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
LINK[k] 0 2 1 12 5 4 11 9 8 10 7 6 14 13 15

で済むのである. 1のLINKが空いているのは, もったいないようにも思うが, 1はこの二進木全体のもっとも左のノードなので, 二進木がスーパーノード0の右部分木であるとすれば, 0の右リンク3を置かせてあげればよく, 従って1のLINKは3とするのが正解である.

k 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
LINK[k] 3 2 1 12 5 4 11 9 8 10 7 6 14 13 15


上の二進木を中順序で辿る通常のアルゴリズムは,

(define (traverse k)
(if (not (eq? (llink k) '())) (traverse (llink k)))
(display k) (display " ")
(if (not (eq? (rlink k) '())) (traverse (rlink k))))

(define (llink k)
(list-ref '(() () 1 () () 4 () () 8 () 7 6 () 13 ())
(- k 1)))
(define (rlink k)
(list-ref '(2 () 12 5 () 11 9 () 10 () () 14 () 15 ())
(- k 1)))

(traverse 3) => 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

リンクを

(define (link k)
(list-ref '(3 2 1 12 5 4 11 9 8 10 7 6 14 13 15) (- k 1)))

とした場合のllinkとrlinkは

(define (llink k) (let ((l (link k)))
(if (< l k) l '())))
(define (rlink k) (let ((l (link (+ (modulo k 15) 1))))
(if (> l k) l '())))

(traverse (link 1)) => 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

なるほど.

2009年11月4日水曜日

ドーナッツを截る

前のblogでドーナッツの断面の平面図を描いたが, ドーナッツがどう見えるかは, あれからでは, よく分からない人が多いのではないか.

そこでもう少し違う角度の絵を描いてみることにした.



この図は, 前の最後の絵に似ている. ただ楕円が2つあり, それぞれ直かに楕円として描いた. 青の楕円は, 長軸の上が(0,0.5), 下が(0,-1.5), 短軸の右が(√3/2,0), 左が(-√3/2,0)なので, 座標の原点を0.-0.5に移動し, スケールを(√3/2,1)に設定して, 半径1の円を描いている.

さて, 任意の鉛直面P,P'でドーナッツを切ったとき, その面と交わる楕円の位置は(y+0.5)2+4x2/3=1とy=kxの交点だから, EとFは簡単にえられる.

P,P'がドーナッツを切る2つの円を母線ということにすると, EやFがその母線と交わる点も分かり, 母線の円の中心からの角度も分かる.

例えばAでは青楕円も緑楕円も120度になり, その点でのドーナッツは切断面と接するので, 母線全体は残る. y'の方向では, 青とは0度, 緑とは180度で, ドーナッツの下半分が残り, 上半分は切り取られる.

これを細かい間隔の母線で決定し, それを描いたのが, 次の図である.



最初の図は, 切り取る部分がまだ置いてあって, 切り取る部分の母線には, 色をつけた. 色はxからy'が赤, y'からx'が青, x'からyが緑, yからxが橙になっている.



もう一つの図が上半分をどけた図で, 切断面は破線で描いてある. 隠面消去していないので, 描いてはみたが, やはりちらついて, 納得出来ていない. 母線の一部の円弧は, みな下向きに垂れ下がっていると思って見る必要がある.

2009年11月1日日曜日

ドーナッツを截る

The New Margin Gardner Mathematical LibraryにSliced Doughnutという話題があった. ドーナッツを平面で截ると切口はどうなるかというのだ. 中心を通る水平面で截ると, 同心円になる. 中心を通る鉛直面で截ると, あなの幅だけ離れた2つの円になる. さらにある面で截ると, 切口は交差する2つの真円になる. これはどういう場合かという問題である. ドーナッツの直径は3インチ, 穴の直径は1インチとする.

立体の切口を想像するのは, 結構難しい. 円錐を截って楕円が現れることなど, 証明がなければ直ぐには信じられない.

まぁこんなことだろうと, 解答を見たら, 一応合ってはいたが, これもあてずっぽだったで, もう少し確かめたいと思った.

立体の断面の作図は, 駒場の学生のころ, 得意中の得意だったので, それをやってみたら, 存外簡単であった.

上の図はドーナッツの平面図と側面図である. 側面図に鉛直面で截ったときの, 断面の2つの円が見えるが, それに接する, 赤線の平面で截ると, 断面が真円になるらしい.



中の図が, 断面の描き方である. 側面図の中央に座標原点があるとする. 高さtの水平面で截ったときの, ドーナッツとの境界が, C0, C3と, C1, C2で, それを半径をする同心円を, 平面図に示した. 高さtの水平面と斜めの断面との交点をxとすると, ドーナッツの切口は平面図のy0, y1, y2やy3を通る.



従って, tをドーナッツの最下面から最上面まで, 少しずつ動かしながら, yの点を求めてつなげればよい. こういう作業にはPostScriptは適している.

下の図がそうやって描いた断面で, 楕円をしている. もう1本は, 上下対称のはずなので, 省略した. A点は, 截る面がドーナッツに接するところで, 平面図においてAから上を通りBまでは, y>0の点をつなぐ. Bを過ぎ, C,D,Aはy<0の点をとる. Cから下を通ってCまでは, 外側の同心円との交点である.




この楕円の上下の長さは2インチである. 赤線の勾配は30度, 赤線を斜辺をする直角三角形の高さは1なので, 斜辺の長さはちょうどは2インチになり, 真円であることが分かる.

2009年10月28日水曜日

ピタゴラ装置

NHKの子供の番組にピタゴラスイッチというのがある. そのピタゴラ装置で因果の列が次々と起きるのがまぁ面白い. 私はプログラムの面白さと相通じるものがあると考えている.

ところで, 11月28日の日経新聞の最後の面に, こういう写真があった.



記事は「転がれコロコロ遊び演出」というので, ピタゴラの仲間のようなものを作る話であった. そして写真のキャプションは「金属の玉が上がっていくようにみえる「コロコロ」」というのである. 新聞に説明はなかったが,なるほど, そういうことかと分かったので, Processingでシミュレートするプログラムを書いた.

http://playground.iijlab.net/~ew/korokoro/korokoro.html

にあるので, 見てください. 画面の内部でマウスをクリックすると, 最初に戻り, 再度始まる. (立ち上がるのに時間が掛るかもしれない.)

そういわれてみれば, 玉が上がっていくようにも見えなくはない.

2009年10月18日日曜日

日光御成道

川口から幸手まで歩いたのに続き, 川口の対岸の岩渕から本郷追分まで歩いた. 川口側は鎌倉橋という御成街道の出発点があるが, 岩渕側はどこが渡河点か分からないので, 新荒川大橋の, 荒川と新河岸川の中間の土手の上から出発した.

それから後は, 北本通り(きたほんどおり)を, 王子, 駒込, 本郷と歩くだけである. この道は自分でも車を運転したり, バスで通ったりしているので, 記録のための消化試合のようなものだ.

志茂のあたりで, 赤羽の辺から東の方に見える, Coit Towerのような塔が突然近くに現れた.

Coit Towerは, サンフランシスコへいった人は覚えていると思うが, テレグラフヒルに立つ塔である. そのCoit Towerにそっくりな塔は何かと思ったら, 北清掃工場であった.

王子から飛鳥山を過ぎたあたりに, 西が丘の一里塚が残っている. 御成街道で江戸から2里目である. 1里目は本郷追分にあったらしいが, 今はなにも残っていないから, これが最初の一里塚だ. (中山道の方は, 志村の一里塚が3里目で残っている.)



西が丘の一里塚


東大農学部の前を過ぎ, 本郷弥生の交差点を過ぎたところに, 日本橋へ4Kmの距離標があるので, そのあたりに一里塚があったのでないか.

東京大学が本郷に終結し始めたのは, 医学部の明治10年(1877年)だ. その頃の地図を見ると, 一里塚があったかもしれない. 夏目漱石の三四郎は, 明治40年頃の話だが, 兼安は出てきても, 一里塚の記録はない.

ところで, 岩渕から追分まで, 歩行2時間であった.

2009年10月13日火曜日

日光御成道

暑くもなく寒くもない街道歩きの季節だ. 川口から幸手まで日光御成道(にっこうおなりどう または日光御成街道)を2回に分けて歩いた. 37Km

そもそも日光御成道は, 徳川の将軍が家康をまつる日光東照宮へ参詣する際, 最初から日光街道(日光道中)を行くのではなく, 途中岩槻城に宿泊するため, 江戸から幸手まで岩槻経由の街道を整備したものである.

日光街道は日本橋から北へ向けて, 中山道と連れ立って出発するが, 室町で右折して分かれてしまう. 中山道は東大農学部の前の本郷追分で左折し, 一方, 御成道はそのまま直進し駒込へ向う. 王子, 志茂, 赤羽を経て岩渕で荒川を渡る. 東大前駅からは, 東京メトロ南北線沿いである.

荒川を渡ると川口. そこから鳩ヶ谷, 大門, 岩槻, 幸手の各宿と進む. 大門あたりまでは, ほぼ埼玉高速鉄道に沿っている.

今回の街道歩き. 一日目は昼から出発. 川口でまず善光寺へ寄ってみた. ここは善光寺号の機関車を陸揚げしたところで, それゆえ機関車は善光寺号と呼ばれることになる. 善光寺は東北線の荒川鉄橋から, 下流左岸の間近に見える青銅の屋根の寺院である.

その先の鎌倉橋の碑のところから, 13時半に街道歩きは始まった. 川口を通過し鳩ヶ谷に向う. 鳩ヶ谷変電所のところから105号線に分岐したので, 交通量は少し減る. 武蔵野線を東川口のあたりで越え, 大門へ. 大門は川口, 鳩ヶ谷に続くこの街道の宿場だが, もっとも宿場らしい落ち着いた通りであった.

やがて東北自動車道をオーバークロス. 今回はこの辺から, 足が少しずつ痛くなり, 野田の鷺山を過ぎた上野田で16時頃, リタイアした. 本来は岩槻までいくはずであった. その交差点から大宮駅東口行きのバスで帰宅した.

一週間後, バスで上野田へ戻り, そこから街道歩きを再開. 前回のように途中で足が痛くなったらどうするか, 心配した. 途中東武野田線の岩槻, 東武伊勢崎線の和戸, 東武日光線の幸手を通るが, それ以外は東武バスか旭バスがわずかに走っているだけだ.

まぁその時はその時とばかり, 歩き続けた. 幸いなことに, 今回は最後まで足は持った. 上野田を10時に出発. 11時半頃岩槻を通過. 65号線を歩いて, 幸手着15時であった.

こういう道でも, 一里塚はいくつか残っていた. 江戸から幸手まで12里なので, 全部あっても12個だが, その一つは一里塚を壊して岩槻区役所(もと市役所)が作られている. 相の原に一里塚が残っているそうだが, 気づかなくて残念であった. 一番ちゃんとしていたのは, 白岡下野田のであった.



膝子の一里塚




白岡下野田の一里塚


埼玉教育委員会編集 歴史の道調査報告書第二集 日光御成道 1984

2009年10月11日日曜日

正方化長方形

Tne New Martin Gardner Mathematical Libraryの2冊目にSquaring the Squareという話題があった. 正方化正方形とでもいおうか. 正方形を相異る正方形で埋め尽す問題である. 同じ正方形で埋め尽すなら22が4とか33が9とかに分割すれば出来るが, 相異るとなると出来るかどうか不明である. 本書によると1936年頃, ケンブリッジ大学の4名の学生が挑戦したらしい.

周囲が正方形なら非常に困難だが, 相異る正方形で長方形を敷き詰めるのは, 比較的容易らしい. 下の図の左がその一例で, 横61, 縦69の長方形が相異る正方形で敷き詰めてある. 正方形の中の数が, 正方形の辺の長さで, 最小のは2と書いてある. こういう図形を正方化長方形(squared rectangle)という.

例えばこの上に61掛ける61の正方形, 横に69掛ける69の正方形を置いても, 正方化長方形になるが, 内部に正方化長方形を持たないものだけを単純正方化長方形ということにする.



目から鱗なのは, これが右のような有向グラフ(Smith Diagram)に変換できることだ. つまり左の図の横線を節点にし, 正方形を弧として描き直すと右の図になる. 弧の値は正方形の辺の長さだが, このグラフの電気回路とみると, 弧の値はこの回路の電流となる. また各線の抵抗を1オームとすると, 弧の値は両端の電圧になる.

このグラフでは, 節点の入力弧は横線の上の正方形の幅の和だし, 出力弧は下の正方形の幅の和なので, 電気回路とすると, 入る電流と出る電流は等しことになる. また例えばAとDの間の電圧はAD直接でも, ABCDと回っても同じである. またこの回路のAへ外から入る電流, Fから外へ出る電流を1とすると, 次のような関係が得られる.



の左の図のように,

電流の出入り

A a + b = 1
B c + d = b
C e + f = c
D g + h = a + e
E i = d + f + h
F 1 = g + i

電圧関係

AD a = b + c + e
BE c + f = d
CE e + h = f
DF g = h + i

従って次の係数の連立方程式を解く.

a b c d e f g h i
A 1 1 1
B -1 1 1 0
C -1 1 1 0
D -1 -1 1 1 0
E -1 -1 -1 1 0
F -1 -1 -1
AD 1 -1 -1 -1 0
BE 1 -1 1 0
CE 1 -1 1 0
DF 1 -1 -1 0

(define eq '(
( 1 1 0 0 0 0 0 0 0 1)
( 0 -1 1 1 0 0 0 0 0 0)
( 0 0 -1 0 1 1 0 0 0 0)
(-1 0 0 0 -1 0 1 1 0 0)
( 0 0 0 -1 0 -1 0 -1 1 0)
( 0 0 0 0 0 0 -1 0 -1 -1)
( 1 -1 -1 0 -1 0 0 0 0 0)
( 0 0 1 -1 0 1 0 0 0 0)
( 0 0 0 0 1 -1 0 1 0 0)
( 0 0 0 0 0 0 1 -1 -1 0)
))

とし, 掃出し法のプログラムを書いて解いてみた. Schemeは分数で計算してくれるので, 最後は(-61/37 -28/37)が残る. つまり-61/37i=-28/37 だからi=(/ -28/37 -61/37) -> 28/61 従ってFから出る電流が1なので, g=33/61となる. 従って上の右の図のように, 33と28が得られる.

右の図はもう少し複雑な回路である.

a b c d e f g h i j k l m
A 1 1 1 1
B -1 1 1 0
C -1 1 1 0
D -1 -1 1 0
E -1 -1 -1 -1 -1
F -1 1 1 -1 0
G -1 -1 1 0
H -1 1 1 0
AD 1 -1 1 -1 0
AF 1 -1 1 -1 0
BE 1 -1 -1 -1 0
CG 1 -1 1 -1 0
FE 1 -1 1 0
HE 1 1 -1 0

(define eq '(
( 1 1 1 0 0 0 0 0 0 0 0 0 0 1)
(-1 0 0 1 1 0 0 0 0 0 0 0 0 0)
( 0 -1 0 0 0 1 1 0 0 0 0 0 0 0)
( 0 0 0 0 -1 -1 0 1 0 0 0 0 0 0)
( 0 0 0 -1 0 0 0 0 0 -1 -1 0 -1 -1)
( 0 0 0 0 0 0 -1 0 1 1 0 -1 0 0)
( 0 0 0 0 0 0 0 -1 -1 0 1 0 0 0)
( 0 0 -1 0 0 0 0 0 0 0 0 1 1 0)
( 1 -1 0 0 1 -1 0 0 0 0 0 0 0 0)
( 0 1 -1 0 0 0 1 0 0 0 0 -1 0 0)
( 0 0 0 1 -1 0 0 -1 0 0 -1 0 0 0)
( 0 0 0 0 0 1 -1 1 -1 0 0 0 0 0)
( 0 0 0 0 0 0 0 0 1 -1 1 0 0 0)
( 0 0 0 0 0 0 0 0 0 1 0 1 -1 0)
))

掃出した最後が(-336/167 -99/167)となり, (i j k l m)が後から順に(5/112 3/14 19/112 9/112 33/112) と得られる. つまり横幅は112の正方化長方形になったわけだ.

(ウェブでこういうページをみつけた.)

2009年9月10日木曜日

個人用電卓のプログラミング

今年の3月頃にこのブログに投稿した, 電卓HHCのプログラミングの話題の続きである. この電卓のシミュレータには, スタックが2個あり, その間でデータが受け渡せる. このスタックを, 両方へ延びるテープ, 命令UとNをヘッドの右, 左への移動と思うと, 2スタックモデルはTuring Machineに非常によく似てくる.

そこで今回は, Turing Machineをシミュレートしてみた.

Turingの1936年の有名な論文の始めの方に, テープに0010110111011110...を書く例がある.

つまり, 0の列の間に1を0個, 1個, 2個, 3個, ...と挟むのである. この論文のテープは, 計算結果を書くますと, 作業用の記号を書くますが, 交互に並んでいる. 従って, 上の出力も,テープ全体で見ると, bを空白のますとして, 0b0b1b0b1b1b0b1b1b1b0b1b1b1b1b0...と出来上がる.

Turingのプログラムは以下のようである. (状態bにPe, つまりeを書くとあるが, Turingの論文では, eを上下反対にした形である.)

state symbol operations next state

b Pe,R,Pe,R,P0,R,R,P0,L,L o

o 1 R,Px,L,L,L o
0 q

q 0or1 R,R q
none P1,L p

p x E,R q
e R f
none L,L p

f 0or1 R,R f
none P0,L,L o

開始の状態はbで, 「eを書き右へ行きeを書き右へ行き0を書き右へ行き右へ行き0を書き左へ行き左へ行き」だから, テープは「ee0 0」となり, ヘッドは左の0のところにあって, 新しい状態はo.

状態oは, 結果のますで, 左へ1が続く限り, その右に作業用のxを書く. 0であったら, 0の場所にいるままで, 状態qへ. 従って, 最初はすぐqになる.

状態qでは, 結果のますに0が1がある限り右へ行き, 空白を見つけたら1を書き, 左へ進み(つまり作業ますに移り), pになる.

pでは, 作業用のますを左に見て生き, xに出会えば, それを消してqになるから, 右の空白に1を書きpに戻る. つまりxの数だけ, 右に1を書く. その前に, oからqに移ってきた時にも1を書いたから, xの数より1個多く1を書く. xが見つからず, 左端のeが見つかれば, 右に行き, 結果のますに移ってfになる.

状態fでは結果のますを右に見ていき, 空白を見つけたら0を書き, 左へ2歩行き, 1の真下にいて, oに戻る.

oは前に述べた通り, 1が続くとxを併記していき, 次回1を書き続ける時の準備をする.

とりあえずSchemによるシミュレーションはを下に示す. テープ長は50で, それを越えて読み書きしようとすれば, プログラムを修了する. 記号e, x, 空白を, 2, 3, 4で表現する.

(define e 2) (define x 3) (define n 4)
(define tapelength 50)
(define tape (make-vector tapelength n))
(define i 0)

(call-with-current-continuation
(lambda (throw)

(define (put x) (if (< i tapelength)
(vector-set! tape i x) (throw 'ok)))
(define (get) (if (< i tapelength) (vector-ref tape i)
(throw 'ok)))
(define (r) (set! i (+ i 1)))
(define (l) (set! i (- i 1)))
(define (eraze) (put n))
(define (b)
(put e)(r)(put e)(r)(put 0)(r)(r)(put 0)(l)(l)(o))

(define (o) ;状態o
(let ((c (get)))
(if (= c 1) (begin (r)(put x)(l)(l)(l)(o))
(if (= c 0) (q)))))

(define (q) ;状態q
(let ((c (get)))
(if (or (= c 0) (= c 1)) (begin(r)(r)(q))
(begin (put 1)(l)(p)))))

(define (p) ;状態p
(let ((c (get)))
(if (= c x) (begin (eraze)(r)(q))
(if (= c e) (begin (r)(f))
(begin (l)(l)(p))))))

(define (f) 状態f
(let ((c (get)))
(if (or (= c 0) (= c 1)) (begin (r)(r)(f))
(begin (put 0)(l)(l)(o)))))

(b)))

(define (printtape) ;テープを文字列にする. (0 0 1 4 2 3)
(let ((t (vector->list tape))) ;=>"001 ex"
(display (list->string
(map (lambda (c)
(cond ((= c 4) #\space) ((= c 0) #\0)
((= c 1) #\1) ((= c 2) #\e)
((= c 1) #\1) ((= c 2) #\e)
((= c 3) #\x))) t)))))

最後のprinttapeは, 記号のリストを文字列に変換する. throwで飛び出した時のテープは

(printtape)
ee0 0 1 0 1 1 0 1 1 1 0 1 1 1 1 0 1x1x1x1 1 0 1 1

qで3つめの1を書こうとした時, 配列の範囲を飛び出したことが分かる.

さて, これを電卓プログラムに書き直す. 元の状態名には, 電卓プログラムでは使えない文字もあるので, b,o,q,p,fをそれぞれg,h,i,j,kとする. そして出来たプログラムは次のとおり. 最初のtはテープに相当する2つのスタックを用意する. スタックの底を示すため, -1を置く. スキャンするヘッドの下には, 空白を置く. rとlはヘッドを右や左に動かす. つまりrはU, lはNなのだが, 底を突き抜けないようなチェックが必要である. PやEもpやxとして実装した.

(1_"U4)t= ;Turing Machine Initialize
(U"4_G2_"G4E)r= ;Move Right
(N"4_G3_"GU4E)l= ;Move Left
(^)p= ;print a letter ap
(4p)x= ;erase print blank

状態b, o, q, p, fはそれぞれサブルーチンg, h, i, j, kとして, 定義する. セミコロンから右は, コメントである. プログラムの上の0や1やnoneは, その条件でジャンプする; またはジャンプからここに飛び込むというような, コメントである.

(2pr2pr0prr0pllh)g= ;b Pe,R,Pe,R,P0,R,R,P0,L,L ->o
; 1 0
("_5_Gi7_"Gr3plllh)h= ;o 1 R,Px,L,L,L -> o 0 ->q
; none none
("_1+7_Grri4_"G1plj)i= ;q 0 or 1 R,R -> q none P1,L ->p
; 01 e x01 e x
("2-c_G"3-d_G"4-d_Gllj9_"Grk3_"Gxri)j=
;p x E,R ->q e R ->f none L,L->p
; none
("_1+7_Grrk5_"G0pllh)k= ;f 0 or 1 R,R ->f none P0,L,L -> o

1を4個, その右に0を1個書いた直後のテープの状態(つまりスタック)は下のとおり.


(-1 2 2 0 4 0 4 1 4 0 4 1 4 1 4 0 4 1 4 1 4 1 4 0
4 1 4 1 4 1 4)(1 4 0 -1)



両端の-1はテープの終りだ. 左スタックの底には最初に書いたeeに相当する22が見える. その後, 空白を示す4を飛ばして読めば, 0010110111011110になっている. 0を書いてからヘッドが左へ2回動いたので, ヘッドの位置はスタックトップの1のところである.

分かってみると, Turingの元のプログラムを, 無駄がないのがよく分かった. また個人用電卓も, 2スタックモデルにすると, こんなことまで出来るということが見えた.

2009年9月7日月曜日

Dudeneyのパズル

Tne New Martin Gardner Mathematical Libraryの2冊目を読んでいたら, Dudeneyのパズルというのに出会った.

正三角形を4つのピースに切り分け, 再配置して正方形を作れ. 難問である. この本は同時に解答を示しているが, それでよい理由はそう自明ではない.

まず解答は次の通り. 正三角形ABCの辺ABに中点D, BCに中点Eをとる. AEをEの方向にEBの長さだけ延長しFとする. AFの中点Gを中心に半径GFの円を描き, CBのBの方向への延長上との交点をHとする. Eを中心に半径EHの円を描き, ACとの交点をJとする. AC上にBEに等しくJKをとる. DとKからJEに下ろした垂線の足をそれぞれLとMとする.

正三角形ABCをJE, DL, KMで切り分け, 再配置すると, 同面積の正方形が得られる.



元の正三角形の1辺の長さを2とする. 三角形の高さは√3なので, 正三角形の面積も√3である. 従って同面積の正方形の1辺の長さは√√3である.

上の作り方で出来た正方形の1辺は√√3であろうか.

正三角形の1辺の長さを2とし, 作り方の図で考えると, AE=√3, AF=√3+1, ∴AG=GF=(√3+1)/2, GE=√3-(√3+1)/2=(√3-1)/2. GEHにおいて, EH=√(GH^2-GE^2) =√((√3+1)^2/4)-(√3-1)^2/4)) =√((3+2√3+1)/4-(3-2√3+1)/4)=√√3=EJ. これが作るべき正方形の1辺である.

Eから辺ACに下ろした垂線の足をNとする. EN=(√3)/2. 三角形JENと三角形JKMは相似で斜辺から比は√√3:1. KM=√3/2/√√3=(√√3)/2. つまり正方形の1辺の半分である.

ところで三角形EDLも, DE=JK(=1)だから, JKMと合同である. DL=(√√3)/2.

従って正方形の上と下の辺はLD+LD=MK+MK=√√3.

左の縦の辺は, LE=MJだから, LE+EM=MJ+EM=EJ=√√3. 右の縦の辺も同様. 従って正三角形と同面積の正方形が得られた.

なるほど.
ところで
√√3=1.31607,
A=(0,0), B=(1,√3), C=(2,0), D=(1/2,√3/2), E(3/2,√3/2), JN=√(√3-3/4), AJ=AN-JNなので J=(3/2-JN,0), K=(5/2-JN,0), ∠EJN=atan(EN,JK)=41.15゜
である. NとKはほとんど同じ位置だが, JN=0.99098ゆえ, Kの方がわずかに右にある.

2009年9月2日水曜日

中西式置換プログラム

前回のブログ(2009年8月27日中西流置換プログラム)の最後に書いたTAOCPのアルゴリズム7.2.1.2-Lは, Lというだけあって辞書式順(lexicographicalorder)で, n項の置換, 例えば(1 2 2 3)の置換

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

を作る. 辞書式順とは, 配列の左端に小数点があるとして, 数値順にソートするものである.

(0 1 2 3)の中西ソートの結果

(perm '(0 1 2 3)) =>
((0 1 2 3) (1 0 2 3) (1 2 0 3) (1 2 3 0) (0 2 1 3) (2 0 1 3)
(2 1 0 3) (2 1 3 0) (0 2 3 1) (2 0 3 1) (2 3 0 1) (2 3 1 0)
(0 1 3 2) (1 0 3 2) (1 3 0 2) (1 3 2 0) (0 3 1 2) (3 0 1 2)
(3 1 0 2) (3 1 2 0) (0 3 2 1) (3 0 2 1) (3 2 0 1) (3 2 1 0))



(define (dec s)
(string->number (string-append "."
(apply string-append (map number->string s)))))

で小数にする.

(dec '(0 1 2 3)) => .0123


先ほどの置換を小数にしてソートすると

(sort (map dec (perm '(0 1 2 3))) <) =>
(.0123 .0132 .0213 .0231 .0312 .0321 .1023 .1032 .1203 .123
.1302 .132 .2013 .2031 .2103 .213 .2301 .231 .3012 .3021
.3102 .312 .3201 .321)


アルゴリズムLの結果の辞書式順
 
                                                                    
(algorithm7212l '(0 1 2 3))=>
((0 1 2 3) (0 1 3 2) (0 2 1 3) (0 2 3 1) (0 3 1 2) (0 3 2 1)
(1 0 2 3) (1 0 3 2) (1 2 0 3) (1 2 3 0) (1 3 0 2) (1 3 2 0)
(2 0 1 3) (2 0 3 1) (2 1 0 3) (2 1 3 0) (2 3 0 1) (2 3 1 0)
(3 0 1 2) (3 0 2 1) (3 1 0 2) (3 1 2 0) (3 2 0 1) (3 2 1 0))

になる. いまはどれも4桁なので, 小数ではなく, 整数と思ってソートするのでもよい.

私のブログの2008年10月28日のに, 切頭八面体の絵があり, {0,1,2,3}の全順列を, 隣り同士の交換で実現したと書いてある. 隣り同士の交換で, 全順列を作るのを, 英語ではPlain Changeという.

下の図がその1例である. 左の24段が, {0,1,2,3}の置換で, x印のところが隣り同士交換した場所である. 一番下のx印の交換をすると, 一番上に戻る.

右上の6段は, 左で0を挿入される{1,2,3}の置換である. これも同じように出来ている.

赤線は0の移動を示す.




x印が連続するように, 中西アルゴリズムを少し修正すると,plain changeのアルゴリズムになる. 前回のプログラムのappendを1つおきに逆にappendするようにすればよい. それが下のrevap (reverse appendのつもり)

(revap '(0 1) '(2 3) '(4 5) '(6 7)) => (0 1 3 2 4 5 7 6)

これを使うとplain changeが出来る.

(define (plainchange ls)
(define (revap . ls)
(cond ((null? ls) '())
((null? (cdr ls)) (car ls))
(else (append (car ls) (reverse (cadr ls))
(apply revap (cddr ls))))))
(define (ins h tl)
(if (null? tl) (list (list h))
(cons (cons h tl)
(map (lambda (l) (cons (car tl) l)) (ins h (cdr tl))))))
(if (null? ls) (list '())
(apply revap (map (lambda (l) (ins (car ls) l))
(plainchange (cdr ls))))))

(plainchange '(0 1 2 3))の結果を切頭八面体のHamiltonian閉路にしたのが



である.

2009年8月27日木曜日

中西式置換プログラム

置換(順列)を計算するプログラムにもいろいろな流儀がある. あるとき中西流(と私が勝手に呼んでいる)というのを知った.

(2000年に身罷った中西正和君について, 私が以前, 「数式処理」に寄稿したものを, googleでサーチしていて, 偶然見つけた.)

中西君のプログラムは, むかし風にM式で書いてあり,

perm[x]
=[null [cdr [x]] -> list [x] ;
t -> mapcon [perm [cdr [x]];
quote [lambda [[y]; insert [nil; car [x]; car [y]]]]
]]
;
insert [x; a; y]
= [null [y] -> list [append [x; cons[a;y]]] ;
t -> cons [append [x; cons [a; y]];
insert [nconc [x; list [car [y]]]; a; cdr [y]]
]]

いま風に書けば (近年はMS変換の出来る人も少ないに違いない.)

(define (perm x)
(if (null? (cdr x)) (list x)
(mapcon (perm (cdr x))
(lambda (y) (insert '() (car x) (car y))) )))

(define (insert x a y)
(if (null? y) (list (append x (cons a y)))
(cons (append x (cons a y))
(insert (append x (list (car y))) a (cdr y)))))

mapconも定義する.

(define (mapcon x f)
(if (null? x) '()
(append (f x) (mapcon (cdr x) f))))

insertは以下のような関数である.

(insert '() 'a '(b c)) => ((a b c) (b a c) (b c a))

これを下請けにして, perm

(perm '(a b c)) =>
((a b c) (b a c) (b c a) (a c b) (c a b) (c b a))


これでもいいのだが, mapconnconcが気に入らず, わたし風に中西流を書いてみた.

(define (perm ls)
(define (ins h tl)
(if (null? tl) (list (list h))
(cons (cons h tl)
(map (lambda (l) (cons (car tl) l)) (ins h (cdr tl))))))
(if (null? ls) (list '())
(apply append (map (lambda (l) (ins (car ls) l))
(perm (cdr ls))))))

入れ子で定義してあるinsを見てみる. 働きは中西insertと同じである.

(define (ins h tl)
(if (null? tl) (list (list h))
(cons (cons h tl)
(map (lambda (l) (cons (car tl) l)) (ins h (cdr tl))))))

(ins 'a '(b c d)) =>
((a b c d) (b a c d) (b c a d) (b c d a))

3行目の(cons h tl)は, 完成したリストの先頭を作る, つまり(a b c d)を作って, 後から出来てきたリスト((b a c d) (b c a d) (b c d a))consする.

後のリストの作り方は, (ins h (cdr tl))が, (cdr tl)つまり(c d)に, hつまりainsしたもの, ((a c d) (c a d) (c d a))のそれぞれに, (car lt)つまりbconsして作る. mapを使う.

tlnilだったら, ((a))を返す. (list (list h))

(define (perm ls)
(if (null? ls) (list '())
(apply append (map (lambda (l) (ins (car ls) l))
(perm (cdr ls))))))

lsnilの場合は後回しにして, ls(a b c d)だったとする. 最後の方の(perm (cdr ls))により, (b c d)perm, ((b c d) (c b d) (c d b) (b d c) (d b c) (d c b))が出来るわけで, このそれぞれにainsし, それを最後にappendする.

lsが段々短くなり, (d)だけになったときを考える. insでつけるのはd, つけられるリストは()で, そのときは, ()perm(())が返り, ((d))が出来る.

それにcinsするから, ((c d) (d c))が出来る.

このようにしてpermが作れた.

ついでだが, TAOCPのアルゴリズム7.2.1.2-LのScheme版は以下の通り.

(define (algorithm7212l as) ;TAOCP V4F2 p.39
(define (interchange ls a b)
(define (list-set! ls n a)
(set-car! (list-tail ls n) a))
(let ((t (list-ref ls a)))
(list-set! ls a (list-ref ls b))
(list-set! ls b t)))
(let ((n 3) (j 0) (k 0) (l 0) (ps '()))
(define (l1) (set! ps (cons (tree-copy as) ps)) (l2))
(define (l2)
(set! j (do ((j (- n 1) (- j 1)))
((or (< j 0) (< (list-ref as j) (list-ref as (+ j 1))))
j)))
(if (>= j 0)(l3) (reverse ps)))
(define (l3)
(set! l (do ((l n (- l 1)))
((< (list-ref as j) (list-ref as l)) l)))
(interchange as j l) (l4))
(define (l4)
(define (loop)
(if (< k l)
(begin (interchange as k l) (set! k (+ k 1))
(set! l (- l 1)) (loop))))
(set! k (+ j 1)) (set! l n) (loop) (l1))
(l1)))

2009年8月25日火曜日

Martin Gardner Library

アメリカの雑誌で創刊が古いのは, ほかにもあろうが, National Geographic Magazine と Scientific American が双璧であろう. インターネットで調べると, National Geographic Magazine は1888年, Scientific American は1845年創刊とある.

私もずいぶん長い間, Scienetific Americanの読者であった. 終戦後, 日比谷にCIE図書館というのがあり, アメリカの雑誌が自由に閲覧できた. どの雑誌を借りても, 一種独特のにおいがした. その中にScientific Americanもあった.

その後, 自由に購入出来るようになってからは, いつのころからか, Scientific Americanは定期講読するようになった. そうなったのは, やはりMartin GardnerのMathematical Gamesが面白いからであったのは, 否めない. もちろん他の記事, 例えばTuring Machineや, von Neumannの自己増殖機械や, Grey Walterのカメの話など, まだ記憶に鮮明に残っている.

それにしても, Mathematical Gamesのコラムは, 毎月圧巻であった. いろいろの思い出がある. Life Gameもその一つ. またHofstadterのGoedel, Escher, Bach の存在も, そのコラムで知った. そのコラムにはGEBの本の表紙の絵, 立方体に糸鋸で3方向から, G, E, Bを切り出したものが出ていた.

1979年のIJCAI(人工知能国際会議)は東京で開催され, 私もなにかのセッションの座長を頼まれたので, 初日の懇親会に参加した. みると, その表紙の絵のTシャツを着ている若者がいる. 「この絵知っている」というと, 彼は隣りを指し, 「これが著者だ」といった. つまり, 彼の隣りにDoug Hofstadterがいたのだ. Scientific Americanの書評は読んだというと, Dougは本を送るよといってくれ, 9月頃Goedel Escher Bachを入手した. (その若者はScott Kim. Scottも後にInversionsという文字を変形するアートの本を送ってきた.)

Scientific Americanには, またMathematical Gamesコラムには, 生涯忘れられない記事がいろいろだ.

さて, そのMathematical Gamesが単行本のシリーズになっていると, 最近聞いたので, 既刊の分をさっそくAmazonで購入した.


The New Martin Gardner Mathematical Library
, Cambridge University Press, 15 volumes.

懐かしい話題がいっぱいつまっている. 1956年12月のHexaflexagonの話が最初のようだ. 私はそれも読んだ覚えがあるので, きっとMathematical Gamesは初めから見たのかも知れない. コラムの最後の1年はGardnerとDoug Hofstadterが隔月に執筆し, 翌年からはHofstadterが毎月書くことになって, Martin Gardnerのコラムはなくなったのであった.

これでしばらく, ほかの仕事がとどこおるかもしれない. またプログラムを書く題材が目の前にちらつき始めた.

2009年8月24日月曜日

投票数

投票数の続きである. 前回はかっこの配置の話だったので, 左かっこと右かっこの数が同じ場合であった. つまり左は右に追い越されないが, 最後は左と右は同票を得た.

得票数が違う場合はどうか. 候補者Qがq票, 候補者Pがp票とって, Qがずーっと優勢を保って敗けなかったとすると, 下の図のようになろう.



(0,0)から出発し, 斜線を越えることなく点(q,p)までくるコースの数Cpqが知りたい. 図に矢じるしで示すように, (q,p)には(q-1,p)から来る場合と, (q,p-1)から来る場合とがあるから,

Cpq = Cp q-1 + Cp-1 q

である. 境界条件を考えると C0 0 = 1, C0 q = 1, Cp q(p > q) = 0 だから, 図の格子点に合わせてCの値を書くと

q 0 1 2 3 4 5
p3 5 14 28
2 2 5 9 14
1 1 2 3 4 5
0 1 1 1 1 1 1


空白の場所は0と思い, 上の図の矢じるしの元に相当する左と下の値を足せばこの表は次々と作れる. 前回のn=3の時の値5は, q=3,p=3のところにある. 1711年にこの三角形に気づいたのは, de Moivreだそうだ. de Moivreの定理というのを, 高校のころ習ったので, 懐かしい名前である.

この表の出来方は, なんとなくPascal 三角形を連想させるが, 果たして関係はあるか. 前回と同じ推論をするなら, (0,0)から(q,p)まで格子点を経由していく全コースはp+qCpである. このうち, 例の斜線を越えるものを修正コースに変更するとすれば, 前回は2nCn-1 であったのと同様,p+qCp-1が, 途中でPの票がQの票を上回る場合である. 従って

Cpq = p+qCp - p+qCp-1        (1)

となる. Cが2種類あって申し訳ないが, 添字をみればどちらか分かるはず.

ところで2項定理を思い出してみると,

mCn = m-1Cn + mCn-1        (2)

とう式があった. これを(1)に適用すれば

Cpq = p+qCp - p+qCp-1

={unfolding}

(p+q-1Cp + p+qCp-1) - (p+q-1Cp-1 + p+qCp-1-1)

={項の並べ替え (A+B)-(C+D) -> (A-C)+(B-D)}

(p+q-1Cp - p+q-1Cp-1) + (p+qCp-1 - p+qCp-1-1)

={folding}

Cp q-1 + Cp-1 q.


至極当然であった. (こういう式の変形は, 手書きでやっていると思うと大変そうだが, 実はエディタを使っているので, まったく大したことではない)

2009年8月23日日曜日

投票数

延び延びの衆議院選挙のただ中になった. 新聞などによると, 自公が退潮で, 民主が攻勢らしい. 私は治承4年の世の中を連想している. 傲れる平家が都落ちし, 木曾義仲や源義経が京都に向っている.

ところでTAOCPにballot numberの話がでていた. AとBの候補の得票を順に書いていき, 一方が他方を一度も越えない系列の数である. もともとはn対のかっこのならべ方のことで, n=3とすると, 3対の正しいかっこの配置は,

()()(), ()(()), (())(), (()()), ((()))

の5通りしかない. 左から見ていき, 左かっこより右かっこが多くなることはない. つまり, 右かっこの得票が左かっこの得票を越えることがないというので, ballot numberと関係する. つまりn=3のballot numberは5だ.

一般に2nCn - 2nCn-1である. この式の導き方が面白い.

0, 1, ..., 5から3個のものをとる組合せは, 辞書式順で

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

だけある. ちなみにこのリストはTAOCPのアルゴリズム7216-Lをschemeにして

(define c '(0 1 2 6 0)) (define cs '())
(define (loop j c)
(if (cddr c)
(if (= (+ (car c) 1) (cadr c))
(cons j (loop (+ j 1) (cdr c)))
(cons (+ (car c) 1) (cdr c)))
'()))
(define (comb)
(set! cs (cons (cddr (reverse c)) cs))
(set! c (loop 0 c))
(if (list-tail c 3) (comb) 'ok))

と定義し, (comb) と起動して (reverse cs)で得た.

この(2 1 0)は左から0,1,2番目に左かっこを, 残りに右かっこをおくということである. 従って ((()))に相当する. (3 1 0)(()())だ. x,yとも0,1,2,3の格子を書き, (0,0)をStartとし, 左かっこがあれば右へ, 右かっこなら上へ進む. かっこは左右3個ずつだから, 最後は(3,3)のGoalに至る.

この格子で, (0,0)から(3,3)への対角線を上へ越えなければ, 正しいかっこの対応が得られる.

下の図の左でその要領をしめす. 0,1,2,3の4つのコースが描いてある. 0と1は対角線を越えないのでOKだが, 2と3はだめだ. 2は())((), 3は)))(((である.




上へ越えたコースがあれば, 上へ越えた最初の点((0,1)から(2,3)への斜線に遭遇した点)で, それまでのコースの縦横を交換する. 例えば2は横, 縦, 縦と進んだが, これを縦, 横, 横に進んだことにする. 残りの横, 横, 縦はそのまま進む. つまり点(1,2)で, 点(2,1)にいたことにし, 残りを進むので, Goalは(3,3)ではなく(4,2)になる.

次の図は, 20通りのかっこの図のコースを, 上の規則で描いたものである. コースは太線で示す. それぞれの図の下の, 2 1 0のようなのは, 左かっこの位置, その次に 5 4のように書いてあるのは, 修正コースでの縦に移動した位置である. これを見ると0,1,...,5から2個取る組合せはすべて1回ずつ現れている.



最初の式に戻ると, 2nCnはStartからGoalに至るすべてのコースの数. 2nCn-1はGoal'に至る修正コースの数であり, その差が正しいかっこの配置の数であったのだ.

2009年8月8日土曜日

複素数用計算尺

複素数用の計算尺があると知って, 例によってその絵を書いてみることにした. 複素数 x+iy の対数は実部が (log (sqrt (+ (* x x) (* y y)))), 虚部が(atan (/ y x))なので, Schemeで実験する. SchemeにはComplex型があるので, こういう時は便利だ.

(* 2+i 3+2i) => 4+7i

(define (clog x y)
(list (log (sqrt (+ (* x x) (* y y)))) (atan (/ y x))))

と定義し

(clog 2 1) => (.8047189562170503 .4636476090008061) ;log 2+i
(clog 3 2) => (1.2824746787307684 .5880026035475675) ;log 3+2i
(clog 4 7) => (2.0871936349478184 1.0516502125483738);log 4+7i

(+ .8047189562170503 1.2824746787307684) => 2.087193634947819
(+ .4636476090008061 .5880026035475675) => 1.0516502125483735

たしかに(clog 2 1) + (clog 3 2) = (clog 4 7) であった.

次にとりあえず x=1 にし, y= -10から1おきに10まで変えながら, clogをとると,

                                                                           
(do ((y -10 (+ y 1))) ((> y 10))
(display (clog 1 y)) (newline))

(2.30756025842063 -1.4711276743037347)
(2.2033596236321267 -1.460139105621001)
... 7行省略
(.3465735902799727 -.7853981633974483)
(0 0)
(.3465735902799727 .7853981633974483)
...7行省略
(2.2033596236321267 1.460139105621001)
(2.30756025842063 1.4711276743037347)

この数値を元に, 拡大や移動しながら, 曲線を描いてみた.


ここまで出来ればあとはPostscriptの出番である. PostScriptによるプログラムは以下のようだ.

/xscale 240 def /yscale 200 def
50 250 translate
/re {x x mul y y mul add sqrt log} def
/im {x y atan dup 180 gt {360 sub} if 100 div} def
/x 1 def /y 0 def re xscale mul im yscale mul moveto
xscale 2 mul 0 rlineto stroke
1 1 10{/x exch def
/y x -10 mul def
re xscale mul im yscale mul moveto
x -10 mul 0.1 x 10 mul{/y exch def
re xscale mul im yscale mul lineto} for
stroke} for

1 0 0 setrgbcolor
/y 1 def /x 0 def re xscale mul im yscale mul moveto
xscale 2 mul 0 rlineto stroke
1 1 10{/y exch def
/x y -10 mul def
re xscale mul im yscale mul moveto
y -10 mul 0.1 y 10 mul{/x exch def
re xscale mul im yscale mul lineto} for
stroke} for


基本の部分を曲線群を以下に示す.

アルファベットで示す各点の複素数と, 座標は

A 1 (0.0 0.9)
B +i (0.0 0.0)
C 1+i (0.15051499 0.45)
D 2i (0.30103 0.0)
E 1+2i(0.349485 0.265650511)
F 2+i (0.349485 0.634349465)
G 5i (0.69897 0.0)

であり, (* 1+i 1+i)=2i, (* 1+2i 2+i)=5i なども読み取れる.

計算尺として使うには, 一方を透明, 他方を不透明の紙に, この図を2枚用意し, 1+2i(E) * 2+i(F)を計算するには, 不透明の図のEに透明のAを重ね, 透明の図のFの下の不透明の図を位置を読むのである. 透明の紙は不透明の紙と平行に動かさなければならない.

なお, 詳しいウェブページはhttp://cs.stmarys.ca/~dawson/sliderule.gifhttp://ci.nii.ac.jp/naid/110000218130/enにある.

2009年7月29日水曜日

手回し機械式計算機

6月15日のブログで, 階差による計算法を書いた時, y=x2+x+41を例に使った. もちろんこれはEulerが1777年に見つけた素数の式である. (xが0から39までに対してyは素数になる. x=40の時は, 402+40+40+12=(40+1)2となり, 素数ではない.)

最近入手した, IEEE Annals of the History of Computing (Vol.31, No.2, April-June 2009)に, Prototype Fragments from Babbage's First Difference Engineという記事があった. 1822年頃, Babbageは階差機関を開発していた. 結局これは, 複製が上野の科学博物館などに展示されている, 職人のJoseph Clementの作った加算器の模型だけが出来たようにいわれている.

それでも見栄えは結構いいので, 情報処理学会の論文賞のメダルなどには, この加算器がレリーフになっている.

その後, 階差機関No.2は, Babbage生誕200年記念で作られ, ロンドンの科学博物館と, カリフォルニアのComputer History Museumで展示されている.

ところで, 上述の記事によると, この加算器以外に, 加算器用に作られた部品が, あちこちの博物館に保存されているということである. 真偽いろいろあるらしく, 試作品のものなど, 怪しいものもあるらしい. そういうものを集める博物館も博物館だが, それを調べる人も調べる人だ.

実は, 私の気になったのは, 記事に引用されているBabbageによる解説である. 1822 年7月3 日に, BabbageがSir Humphry Davyに送った手紙の一節が示されており, 「it proceeded to make a table from the formula x2+x+41.」と書いてあるのだ.

Eulerが発見してから, 半世紀も経っているので, Babbageもこの式を知っていたのであろう. やはりこれは特別な式なのだ.

2009年7月25日土曜日

長大語計算

TAOCPにbroadword computingという話題がある. 長大なビット列に対する演算法である.

例えば x=1110111101100111 の中で 0111 というパターンを探す問題である.
その答は q=0001000000001000 である. つまりqの1に対応するxの場所が0111の0であり, その右に111があるのだ.

方法は q=¬x & x<<1 & x<<2 & x<<3 とすればよい. 途中の状況を示すと下のようだ.

x =(1 1 1 0 1 1 1 1 0 1 1 0 0 1 1 1)
¬ x =(0 0 0 1 0 0 0 0 1 0 0 1 1 0 0 0)
x<<1 =(1 1 0 1 1 1 1 0 1 1 0 0 1 1 1 0)
x<<2 =(1 0 1 1 1 1 0 1 1 0 0 1 1 1 0 0)
x<<3 =(0 1 1 1 1 0 1 1 0 0 1 1 1 0 0 0)
q =(0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0)

つまり探すパターンの各ビットが0なら¬x, 1ならxのままを, 左シフトで, 答の位置に集め, 全体の&をとれば得られるわけである.

応用問題として, 5月8日のブログの最初にあるQ0とQ1から11011, つまり双子の素数が並んでいる場所を探してみる.

(define x #x2196820D864A4C32816D129A64B4CB6E)

として, x & x<<1 & (¬x)<<2 & x<<3 & x<<4 で良い筈である.

から

が得られる. つまり, 5と7と11と13, 11と13と17と19で, この辺はすぐ思いつく. その次は101と103と107と109であり, さらに191と193と197と199が続く. これだけ素数が連続するから, 他の小さい素数はどこに隠れるか見てみると

101 = 101
102 = 17 3 2
103 = 103
104 = 13 2 2 2
105 = 7 5 3
106 = 53 2
107 = 107
108 = 3 3 3 2 2
109 = 109

のように, 7, 13, 17がうまく2や5と一緒になる.

191 = 191
192 = 3 2 2 2 2 2 2
193 = 193
194 = 97 2
195 = 13 5 3
196 = 7 7 2 2
197 = 197
198 = 11 3 3 2
199 = 199

も同様.

ついでにQ0からQ7までを使い

として, もう少し先まで探すと qが得られ, 最左の1の位置を求める.
(ここでi2bは整数を0と1のビット列にする関数. iandなどは, 2つの整数のビット毎ANDなど, ilshは整数のビット左シフトである.)


(length (member 1 q)) => 940

従って並ぶ双子の最大の素数は
(- (* 940 2) 1) => 1879

であり, この辺りは

1871 = 1871
1872 = 13 3 3 2 2 2 2
1873 = 1873
1874 = 937 2
1875 = 5 5 5 5 3
1876 = 67 7 2 2
1877 = 1877
1878 = 313 3 2
1879 = 1879

次は

(length (member 1 (list-tail q (- 1024 940 -1)))) => 745
(- (* 745 2) 1)) => 1489

1481 = 1481
1482 = 19 13 3 2
1483 = 1483
1484 = 53 7 2 2
1485 = 11 5 3 3 3
1486 = 743 2
1487 = 1487
1488 = 31 3 2 2 2 2
1489 = 1489

(length (member 1 (list-tail q (- 1024 745 -1)))) => 415
(- (* 415 2) 1) => 829

821 = 821
822 = 137 3 2
823 = 823
824 = 103 2 2 2
825 = 11 5 5 3
826 = 59 7 2
827 = 827
828 = 23 3 3 2 2
829 = 829

(length (member 1 (list-tail q (- 1024 415 -1)))) => 100

となり, 先ほどの199, 197, 193, 191が得られる.

2009年7月2日木曜日

手回し機械式計算機

6月15日のブログで触れたComrieの論文には, Duplaを使った第2階差をとる方法が書いてある. それは次のようだ.

等間隔の引数に対する関数値が, a, b, c, d, ...だったとする.
第1階差は, b-a, c-b, d-c, ...だから,
第2階差は(c-b)-(b-a)=a-2b+c, (d-c)-(c-b)=b-2c+d, ... である.

この種の計算機では, 数値を入力するのが厄介(というより, 操作ミスの原因)なので, 入力回数を極力少なくするのに気を使った. 以下の方法では, 「なにを置く」は, 各値について1回だけだ. 第2階差は「記録」のところで得られる.

S R1 R2
aを置く a 0 0
1回転 a a a
bを置く b a a
-3回転 b a-3b a-3b
R2帰零 b a-3b 0
1回転 b a-2b b
1 cを置く c a-2b b
2 1回転 c a-2b+c b+c
3 記録 c 階差 b+c
4 -4回転 c a-2b-3c b-3c
5 R1帰零 c 0 b-3c
6 1回転 c c b-2c
1 dを置く d c b-2c
...

この「cを置く」から, 「dを置く」の直前までが1サイクルである. ここから, R1とR2の仕事を交換し, R2に次の階差が得られる. 以下同じ.

実に見事な方法だ. 次の入力数を置数レジスタに置いたまま, 2つのレジスタを使い, 次の階差と次の次の階差の準備を進めている. Duplaでは, 一方に足しながら, 他方から引くことも出来るが, そのようにしていないのは, 操作ミスを心配したからと思われる.

これは今の言葉でいえば, プログラムを作ったようなものである. 手順を書き下し, 後はこの通りに手を動かすだけでよい. 最初に書いたように, それぞれの値が1回しか, 入力されないのが, 最大のメリットである.

Leslie John Comrie(1893-1950)は, 化学を専攻したが, 後に天文学者になり, 各種の計算機を使って数表を作る仕事に尽力した. 計算機も単に手回し機にとどまらず, パンチカードを使った会計機等も利用した. 「A Computer Perspective計算機創造の歴史」には, Comrieは何回も登場する.

Comrieの写真はここにある.

2009年7月1日水曜日

再帰曲線

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

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



dragon曲線は紙を半分に折り, 折り目を直角に開くのだが, それの相当するのは, dragon曲線を表わす矢印の先に, 同じ矢印を右から合わせることである. 従って1次のdragon曲線の矢印1は, 0次の矢印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倍になっていたので, 結局同じように増殖していたことが判明した.

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

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に八進表現が出来ている.

2009年5月28日木曜日

pi関数

TAOCPにpi関数というのが登場する. 真理値表が 1100 1001 0000 1111 となる4変数の関数である. これは, 種明かしがしてあるが, 円周率πを二進で表わしたものだ. ランダムな真理値表が欲しいとき, πや素数表を使う.

十進のπの値は良く知られている.

TAOCPの第1巻の終りの方に八進のπの値が書いてあり,

3.11037 55242 10264 30215 14230 63050 56006 70163 2112

これを二進にすれば

11.001 001 000 011 111 101 101 010 100 010 ...

だから, 始めの方はたしかに上の真理値表だ.

ところが512ビットの真理値表というのも現れる.

1100100100011...00101

だというのである. 本当かな?と思って計算してみた. ついでだから, 1024ビットまで求めて, 十六進で書いたのが下のものだ.


2行目の最後が512桁で, 45だから00101だった.

この十六進の文字の分布も計算した. 下のリストの左から0, 1, 2,... , fの出現数で, 合計は256になっている.

(apply + '(15 16 17 16 21 17 21 10 14 13 14 19 18 14 17 14))
=> 256

ところでこの計算法だが, Lispのbignumを活用する. 十進法のπの値はspigotアルゴリズムで計算した. もっともここの段階で十六進に出来るはずだが, 面倒だったので, 十進の値を出した. 1024ビットは十進では310桁くらいなので, 320桁を用意した.


小数点は下から319桁にあるから, 10の319乗を用意する.

(define pow10 (expt 10 319))

とする. また最初を1100にしたいから, piを4倍する.

(set! pi (* pi 4))

十六進の印字プログラムも必要

(define (hexprint n)
(define hex "0123456789abcdef")
(display (string-ref hex n)))

そして計算する.

(hexprint (quotient pi pow10))
(set! pi (modulo pi pow10))

ここで最初のcが出力される.

(do ((i 0 (+ i 1))) ((= i 256))
(hexprint (quotient (* pi 16) pow10))
(set! pi (modulo (* pi 16) pow10)))

ここから90fdaa...が得られた.

折角プログラムがあるのだから, 自然対数の底のeでも同じ計算をしよう.


TAOCPの八進の値は

2.55760 52130 50535 51246 52773 42542 00471 72363 6166

で, 十六進1024ビットは以下のとおりである.