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倍になっていたので, 結局同じように増殖していたことが判明した.

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