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にある.