>>> a = [1,2] >>> b = [1] >>> b += [2] >>> a [1, 2] >>> b [1, 2] >>> a == b True >>> a is b False
==演算子は中身を、is演算子は参照を比べています。Javaだと==で参照の同一性、equalsで同値性だったはずなんで逆ですね。
>>> c = [1] >>> d = c >>> c is d True >>> c = c + [2] >>> c [1, 2] >>> d [1]
c = c + [2]とした場合、新しいオブジェクトが作られ、それにcが結合されます。cがもともと結合されていたオブジェクトは書き換えられません。
>>> e = [1] >>> f = e >>> e is f True >>> e += [2] >>> e [1, 2] >>> e [1, 2]
拡張代入の場合、cが結合されているオブジェクト自体が書き換えられます。
>>> for i in xrange(10000000): ... e += [i] ... >>> e is f True
Pythonのgenerator(2) グラフ探索(1)のつづき。
こんな感じでノードクラスというものを作ってみます。フィールドは自分のノード名と隣接するノードインスタンスのリストとしました。
class Node:
def __init__(self, name):
self.name = name
self.nodelist = []
def add(self, list):
self.nodelist += list
このクラスを使って先日のものと同じグラフを作ると
n1 = Node("a")
n2 = Node("b")
n3 = Node("c")
n4 = Node("d")
n5 = Node("e")
n6 = Node("f")
n1.add([n2, n6])
n2.add([n1, n3, n5])
n3.add([n2, n4, n5, n6])
n4.add([n3])
n5.add([n2, n3, n6])
n6.add([n1, n3, n5])
これに対する探索用ジェネレータはこんな感じ。再帰はなくしました。深さ優先はスタックを、幅優先はキューを使うと簡単。深さ優先のほうで隣接するノードのリストをひっくり返しているのは先日と同じ出力を得るためです。
def graph_depth2(node):
route = []
stack = []
stack.append(node)
while stack:
n = stack.pop()
if n not in route:
for l in n.nodelist[::-1]:
if l not in route:
stack.append(l)
yield n.name
route.append(n)
def graph_breadth2(node):
route = []
queue = []
queue.append(node)
while queue:
n = queue.pop(0)
for l in n.nodelist:
if l not in route and l not in queue:
queue.append(l)
yield n.name
route.append(n)
実際に探索させてみると
>>> for t in graph_depth2(n1):
print t,
a b c d e f
>>> for t in graph_depth2(n3):
print t,
c b a f e d
>>> for t in graph_breadth2(n1):
print t,
a b f c e d
>>> for t in graph_breadth2(n3):
print t,
c b d e f a
先日と同じ結果に。
もう一個のほうでも
n1 = Node(1) n2 = Node(2) n3 = Node(3) n4 = Node(4) n5 = Node(5) n6 = Node(6) n7 = Node(7) n8 = Node(8) n9 = Node(9) n10 = Node(10) n11 = Node(11) n1.add([n2, n5]) n2.add([n1, n3, n5]) n3.add([n2, n4, n9]) n4.add([n3, n5]) n5.add([n1, n2, n4, n6, n7]) n6.add([n5, n9]) n7.add([n5, n8]) n8.add([n7]) n9.add([n3, n6, n10]) n10.add([n9, n11]) n11.add([n10])
>>> for t in graph_depth2(n1):
print t,
1 2 3 4 5 6 9 10 11 7 8
>>> for t in graph_depth2(n9):
print t,
9 3 2 1 5 4 6 7 8 10 11
>>> for t in graph_breadth2(n1):
print t,
1 2 5 3 4 6 7 9 8 10 11
>>> for t in graph_breadth2(n9):
print t,
9 3 6 10 2 4 5 11 1 7 8
当然同じ結果に。
先日は距離を導入してダイクストラでもやってみようと思ったんですが、ほとんど変わらなそうなので止めました。リストに追加する際に始点からの距離を計算して、ポップする時に一番始点から近いノードを拾って来ればいいだけですし。
ジェネレータのメリットというのは簡潔さとメモリ使用量の削減にあるので、再帰を使ったり、複雑な制御構造を持たせたりするのは避けた方がよいんでしょうね。
をちょっと考えてみます。まずいかにもアルゴリズムの教科書に載ってそうなグラフから。グラフの表現はとりあえずディクショナリを使って、
graph = {"a":["b","f"],"b":["a","c","e"],"c":["b","d","e","f"],
"d":["c"],"e":["b","c","f"],"f":["a","c","e"]}
こんな感じで。ノードとそれに接続するノードのリストが要素。リストの左側のものが優先的に探索されます。ディクショナリまるごと渡したらジェネレータの利点がなくなるんじゃ、って感はありますが。
import copy
def graph_depth(node, graph):
g = copy.deepcopy(graph)
def graph_depth_inner(node):
if node in g:
yield node
nodes = g.pop(node)
for n in nodes:
if n in g:
for i in graph_depth_inner(n):
yield i
for n in graph_depth_inner(node):
yield n
def graph_breadth(node, graph):
g = copy.deepcopy(graph)
def graph_breadth_inner(queue):
q = []
for n in queue:
if n in g:
yield n
q += g.pop(n)
if q:
for n in graph_breadth_inner(q):
yield n
for n in graph_breadth_inner([node]):
yield n
graph_depthが深さ優先、graph_breadthが幅優先。木を探索する場合にはノード、左部分木、右部分木とはっきり分割できるのでジェネレータを再帰的に呼び出すだけでいいのですが、グラフの場合はどこを探索したのかを覚えておく必要があります。ディクショナリをコピーしてそれを破壊しながら探索するのが楽。探索済みのノードはディクショナリから削除することで判定するようにします。それぞれの関数の内部で定義された**_innerが実際の探索を行います。
ノードaから深さ優先探索する場合、aを出力、graph-{"a"}をノードbから深さ優先探索、bを出力、graph-{"a","b"}をノードcから深さ優先探索...と再帰的にジェネレータを適用できます(graph-{"a"}はgraphからノードaを取り除いたグラフとします)。
>>> for t in graph_depth("a",graph):
print t,
a b c d e f
>>> for t in graph_depth("c",graph):
print t,
c b a f e d
幅優先の場合、探索を開始したノードから等しい高さ(辺数)にあるノードをひとまとめにして再帰。もっとスマートにできそうですが。
>>> for t in graph_breadth("a",graph):
print t,
a b f c e d
>>> for t in graph_breadth("c",graph):
print t,
c b d e f a
もう少し大きいグラフを
graph2 = {1:[2,5],2:[1,3,5],3:[2,4,9],4:[3,5],
5:[1,2,4,6,7],6:[5,9],7:[5,8],8:[7],
9:[3,6,10],10:[9,11],11:[10]}
深さ優先
>>> for t in graph_depth(1,graph2):
print t,
1 2 3 4 5 6 9 10 11 7 8
>>> for t in graph_depth(9,graph2):
print t,
9 3 2 1 5 4 6 7 8 10 11
幅優先
>>> for t in graph_breadth(1,graph2):
print t,
1 2 5 3 4 6 7 9 8 10 11
>>> for t in graph_breadth(9,graph2):
print t,
9 3 6 10 2 4 5 11 1 7 8
少しだけいじると
def graph_depth(node, graph):
g = copy.deepcopy(graph)
def graph_depth_inner(node):
if node in g:
yield node
nodes = g.pop(node)
for n in nodes:
if n in g:
for i in graph_depth_inner(n):
yield i
yield node
for n in graph_depth_inner(node):
yield n
深さ優先で行って帰ってきたり、
>>> for t in graph_depth(1,graph2):
print t,
1 2 3 4 5 6 9 10 11 10 9 6 5 7 8 7 5 4 3 2 1
>>> for t in graph_depth(9,graph2):
print t,
9 3 2 1 5 4 5 6 5 7 8 7 5 1 2 3 9 10 11 10 9
def graph_breadth(node, graph):
g = copy.deepcopy(graph)
def graph_breadth_inner(queue):
q = []
for n in queue:
if n in g:
yield n
q += g.pop(n)
if q:
yield "|"
for n in graph_breadth_inner(q):
yield n
for n in graph_breadth_inner([node]):
yield n
幅優先で同じ高さ(辺数)にあるノードを区切ったりできます。
>>> for t in graph_breadth(1,graph2):
print t,
1 | 2 5 | 3 4 6 7 | 9 8 | 10 | 11 |
>>> for t in graph_breadth(9,graph2):
print t,
9 | 3 6 10 | 2 4 5 11 | 1 7 | 8 |
次はグラフの各ノードをクラスとして、ノード間の距離も考えてみます。あんまり変わらなそうですが。
ジェネレータというものを試してみようと思った。
通常の関数の場合、
def hello():
return 'hello,world!'
>>> hn = hello()
>>> hn
'hello,world!'
実行すると文字列が返ってきます。その結果、hnは文字列を指すことになります。
ジェネレータだと、
def hellogen():
yield 'h'
yield 'e'
yield 'l'
yield 'l'
yield 'o'
yield ','
yield 'w'
yield 'o'
yield 'r'
yield 'l'
yield 'd'
yield '!'
>>> hg = hellogen()
>>> hg
<generator object at 0x00D2CB20>
なんかオブジェクトが作られたようです。関数定義内にyieldが含まれる場合、それはジェネレータとして扱われます。ジェネレータは実行されるとgenerator objectを生成し、それがhgの参照となっているようです。このgenerator objectに対して、
>>> hg.next() 'h'
とnextメソッドを実行するとhが返ってきます。これは一行目のyieldの引数ですね。
>>> hg.next() 'e' >>> hg.next() 'l' >>> hg.next() 'l' >>> hg.next() 'o'
next()を続けると、そのたびに一文字ずつ、別の文字が返ってきます。これらはhellogen()の中でyieldの引数として与えたものです。generator objectは、next()が実行されるたびにyieldの引数をひとつずつ順番に返します。generrator objectは「次によばれたとき、どこから開始すればよいか」を記憶しているようです。
>>> hg1 = hellogen() >>> hg2 = hellogen() >>> hg1.next() 'h' >>> hg1.next() 'e' >>> hg2.next() 'h' >>> hg1.next() 'l'
hellogen()によってふたつのgenerator objectを作ってみます。それぞれの「次の開始位置」は別々に記憶されているのがわかります。hg1とhg2は別のインスタンスだということです。
>>> hg1 is hg2 False
上のhellogen()によって作られたgenerator objectは、「12個の文字を順番に返すという約束」を持っています。12個を超えて文字を取り出そうとするとエラーになります。
...
>>> hg.next()
'!'
>>> hg.next()
Traceback (most recent call last):
File "<pyshell#40>", line 1, in <module>
hg.next()
StopIteration
ジェネレータはfor文で回してやることができます。
>>> for t in hellogen(): print t, h e l l o , w o r l d !
これは、
>>> for t in hello(): print t, h e l l o , w o r l d !
と同じように見えます。つーか結果としては同じで、文字が用意されるタイミングが違います。hello()のほうは、forループに入った瞬間に'hello,world!'の文字列全体が用意され、そこから一文字ずつ取り出してprintしています。hellogen()のほうは、まずhellogen()によって'h'一文字だけが用意され、それを取り出して表示、また次の文字が用意され、それを取り出して表示、と一文字ずつ取得と表示を行っています。
極端な例だと
def gentest(n):
i = 0
while i < n:
yield i ** 2
i += 1
def test1():
for t in [x ** 2 for x in range(1000000)]:
if (t % 10000000000 == 0):
print t
return
def test2():
for t in gentest(1000000):
if (t % 10000000000 == 0):
print t
return
こんな感じで大量に調べなければならない値があってかつ候補のほとんどがフィルタされて消えるときには、range()で候補のリストをまとめて用意するのはメモリの無駄ってことでジェネレータが有効ですね。
上のgentest(n)の引数をとっぱらうと
def gentest():
i = 0
while 1:
yield i ** 2
i += 1
0, 1, 4, 9, 16, ...と無限に要素を返し続けるジェネレータが作れます。
定番のFibonacci級数
def fibs():
a = 0
b = 1
while 1:
yield a
a, b = b, a+b
>>> for t in xrange(15):
print f.next(),
0 1 1 2 3 5 8 13 21 34 55 89 144 233 377
Schemeでいうstreamのようなものなので、SICPを真似てちょっと不思議なこともできる。
def fibs2():
yield 0
yield 1
f = fibs2()
g = fibs2()
g.next()
while 1:
yield f.next()+g.next()
>>> f2 = fibs2()
>>> for t in xrange(15):
print f2.next(),
0 1 1 2 3 5 8 13 21 34 55 89 144 233 377
グラフ探索とかもできそう。
一昨日のエントリはひどいものだった。僕はこれをPythonでやりたかったんだ。
foldl(op, v, [a, b, c, d]) = (((v op a) op b) op c) op d
foldr(op, v, [a, b, c, d]) = a op (b op (c op (d op v)))
そしてこう書いた。
def foldl(proc, elem, lis):
lis.insert(0, elem)
return reduce(proc, lis)
def foldr(proc, elem, lis):
lis.append(elem)
lis.reverse()
return reduce(lambda x, y : apply(proc, [y, x]), lis)
fl = foldl(lambda x, y : x / y, 2.0, [1.0, 2.0, 4.0, 8.0])
fr = foldr(lambda x, y : x / y, 2.0, [1.0, 2.0, 4.0, 8.0])
>>> fl
0.03125
>>> fr
0.5
こんなアホみたいなことをする必要なんて全くなかったんだ……。もう少し考えるべきだった、Pythonがこんな不便なはずはないって。
import operator
def foldl(proc, elem, l):
lis = l[:]
lis.insert(0, elem)
return reduce(proc, lis)
def foldr(proc, elem, l):
lis = l[::-1]
lis.insert(0, elem)
return reduce(lambda x, y : proc(y, x), lis)
fl = foldl(operator.truediv, 2, [1, 2, 4, 8])
fr = foldr(operator.truediv, 2, [1, 2, 4, 8])
>>> fl
0.03125
>>> fr
0.5
これでよかったんじゃないか……。そして忘れていた、普通の言語にはループってものがあるってことを。
def foldl_loop(proc, elem, l):
lis = l[:]
for x in lis:
elem = proc(elem, x)
return elem
def foldr_loop1(proc, elem, l):
lis = l[::-1]
for x in lis:
elem = proc(x, elem)
return elem
def foldr_loop2(proc, elem, l):
lis = l[:]
while lis:
elem = proc(lis.pop(),elem)
return elem
時間を計ってみよう。largelist = range(10000000)として、リストの整数を全て足しあわせてみる。
foldl(operator.add, 0, largelist) 3.36605238934 foldl_loop(operator.add, 0, largelist) 5.12172161546 foldr(operator.add, 0, largelist) 7.88937811929 foldr_loop1(operator.add, 0, largelist) 5.67820882263 foldr_loop2(operator.add, 0, largelist) 8.93474310282
reduceはえー。そして気づいた、こんなことしたってHaskellの勉強にはならないって。
foldl,foldr
Hugs.Base> foldl (/) 2 [1, 2, 4, 8] 0.03125 Hugs.Base> foldr (/) 2 [1, 2, 4, 8] 0.5
foldl (/) 2 [1, 2, 4, 8] は、(((2 / 1) / 2) / 4) / 8
foldr (/) 2 [1, 2, 4, 8] は、1 / (2 / (4 / (8 / 2)))ということですね。
Hugs.Base> foldl (/) 3 [1, 3, 5, 7, 2, 1] 0.0142857142857143 Hugs.Base> foldr (/) 3 [1, 3, 5, 7, 2, 1] 1.42857142857143
Schemeで書いてみる。Schemeのライブラリ、SRFIのfoldは
> (fold / 2 '(1 2 4 8)) 8
これは8 / (4 / (2 / (1 / 2)))になってしまう。自分で書くしかない。
(define (foldl proc elem list)
(define (foldl-iter list ans)
(cond ((null? list) ans)
(#t (foldl-iter (cdr list) (proc ans (car list))))))
(foldl-iter list elem))
(define (foldr proc elem list)
(define (foldr-rec list)
(cond ((null? list) elem)
(#t (proc (car list) (foldr-rec (cdr list))))))
(foldr-rec list))
> (foldl / 2.0 '(1 2 4 8))
0.03125
> (foldr / 2.0 '(1 2 4 8))
0.5
> (foldl / 3.0 '(1 3 5 7 2 1 ))
0.014285714285714287
> (foldr / 3.0 '(1 3 5 7 2 1 ))
1.4285714285714284
いい感じに頭がこんがらがってきました。
次は Pythonで。
def foldl(proc, elem, lis):
lis.insert(0, elem)
return reduce(proc, lis)
def foldr(proc, elem, lis):
lis.append(elem)
lis.reverse()
return reduce(lambda x, y : apply(proc, [y, x]), lis)
fl = foldl(lambda x, y : x / y, 2.0, [1.0, 2.0, 4.0, 8.0])
fr = foldr(lambda x, y : x / y, 2.0, [1.0, 2.0, 4.0, 8.0])
>>> fl
0.03125
>>> fr
0.5
>>> print foldl(lambda x, y : x / y, 3.0, [1.0, 3.0, 5.0, 7.0, 2.0, 1.0])
0.0142857142857
>>> print foldr(lambda x, y : x / y, 3.0, [1.0, 3.0, 5.0, 7.0, 2.0, 1.0])
1.42857142857
Pythonて絶対こんな書き方をする言語ではないと思った。
forcedethがnForce410を認識してた
昨日のエントリで、etchにしたらネットに繋がらなくなった、と思ったらやっぱ繋がったと書きました。
でも今朝起動してみたらやっぱり繋がりませんでした。ちょっと調べて原因は判明しました。あまりにバカバカしいことだったんですが書いておきます。
$ dmesg | grep eth forcedeth.c: Reverse Engineered nForce ethernet driver. Version 0.56. forcedeth: using HIGHDMA eth0: forcedeth.c: subsystem: 01509:6006 bound to 0000:00:14.0 eth1: VIA Networking Velocity Family Gigabit Ethernet Adapter eth1: Ethernet Address: 00:02:2A:DD:30:9E eth0: no IPv6 routers present
つまりNICが二つあることが原因だったわけですね。PCのチップセットにはもともとnForce410が組み込まれてたんですが、sargeではこれを認識してくれなかったためGbE-PCI2を新しく買って来てそっちにケーブルを繋いで使ってたわけです。でもetchになってforcedethがnForceを認識するようになり、nForceのほうにeth0が割り振られるようになったと。僕はそんなことも知らずにケーブルの繋がっていないeth0から必死でDHCP DISCOVERを送ったり、固定IPにしてPing撃ってみたりしてたわけです。ケーブルをnForce側に差し替えたらすぐ繋がりました。ひととおり環境が整ったら性能測定でもしてどちらを使うか決めようと思います。
ところでチップセットとカードの場合も二枚差しというのでしょうか。カードを二枚差すのが二枚差し？
Ant
Apache Ant ユーザマニュアル/Apache Ant User Manual
# aptitude install ant
Eclipse
# aptitude install eclipse
C
# aptitude install build-essential
DrScheme
# aptitude install drscheme
Gauche
# aptitude install gauche
Python2.5
# aptitude install python2.5
Haskell
# aptitude install ghc hugs
Ogawa::Memoranda FeedBurnerに移行した件について。を参考にさせていただき、登録&リダイレクトの設定を行う。
新しいフィードはこちら。http://feeds.feedburner.com/panopticon/pGXt
ところでPing-o-Matic!というものを初めて知った。
に申し込んだ。日本語が正しく読めれば（これが難しいのだが）そこそこはいけるはずだ。
参考書は日経が出している完全教本を選んだ。16、17、18年度の過去問が付録としてついているのでやたら分厚くなってる。これ過去問部分だけ取り外せればなー。ちょっと重いので講義のある日には持ち運びたくない。
情報処理技術者試験はいちおう国家資格なので受験料が安くていいのだが、ちょっと勉強しようと思うと受験料以上に参考書代が高くなってしまったりする。それってどうなのって感じはするけど落ちれば受験料は無駄になるし、高度資格の場合は来年まで待たなくてはならない。それらを勘案すると結局アイテックの本も買わないわけにはいかないだろう。僕は参考書に書き込みをするタイプなので売るわけにもいかず、いつしかベッドの脇には参考書類がうずたかく積み上げられてゆくのであった。
Debian GNU/Linux 4.0 (`etch') リリースノート (Intel x86 用) を参考に。
sources.listを書き換える
sources.listの指定をsargeからetchに変更。
deb http://ftp.jp.debian.org/debian/ etch main deb-src http://ftp.jp.debian.org/debian/ etch main deb http://ftp.jp.debian.org/debian/ etch main contrib non-free deb-src http://ftp.jp.debian.org/debian/ etch main contrib non-free
ディストリビューションのupgrade
# aptitude dist-upgrade
けっこう時間かかる。高校のころの参考書が一冊読めた。
GRUBの設定
kernelも新しいのが入るようだったのでやっとく。
# update-grub
Xorgのインストール
パッケージからインストールする場合、ほとんどXFree86と変わらない。Gnomeも入れる。かなり起動が早くなってる。ログインマネージャがどっちかと一緒に入ってた。
# aptitude install x-window-system-core # aptitude install gnome
設定を変更するには、
# dpkg-reconfigure xserver-xorg
日本語環境
# aptitude install kinput2-canna im-switch
localeがいつの間にかLANG=Cになってた。このままだとim-switchが動かないみたいなのでja_JP.EUC_JPに直す。
# im-switch -s kinput2-canna `xinput-ja_JP' を提供するために `/etc/X11/xinit/xinput.d/kinput2-canna' を使います。
これだけで動くようだ。入力メソッドの切替えもこのコマンドで行えるらしい。便利。
フォントの設定
# aptitude install defoma x-ttcidfont-conf # aptitude install ttf-kochi-gothic ttf-kochi-mincho # aptitude install ttf-sazanami-gothic ttf-sazanami-mincho # aptitude install ttf-vlgothic ttf-mikachan
Firefoxのインストール
# aptitude install mozilla-firefox
とするとFirefoxではなくIceweaselというブラウザがインストールされる。名前は違うけど中身はFirefoxと同じものらしい。weaselってイタチか。
かわいすぎ。
Debian Project, Firefoxの名称を変更する方向へ
再起動
してみたらネットにつながらなくなってた。どうやらDHCPまわりがあやしい、と思ってWindowsのほうで解決策を調べる。これでいけそうだ、という手段が見つかり意気陽々と戻って来ると、なんもいじってないのに今度はネットにつながった。あれ？でも設定だけは直しておくか、と思いconfをみてみると、直そうと思ってた部分はどこかの親切なソフトがすでに書き込んでくれてた。畜生。
Java
# aptitude install sun-java5-jdk sun-java5-demo sun-java5-source sun-java5-doc sun-java5-plugin sun-java5-fonts
ドキュメントも入れろや、みたいなメッセージが出るので、Java SE Downloadsへ行ってJ2SE 5.0 Documentationを取得。日本語英語どっちかでよい。sun-java-6-*のほうをインストールした場合は、日本語ドキュメントがないみたいなので英語を入れる。ダウンロードして/tmpに持って行けば、あとはaptが勝手にやってくれる。
ちょっと暇ができたので、scheme分を補充しておこう……。Durand-Kerner法(DK法)というのはニュートン法の改良版で、代数方程式の全ての解を一度に求めることができる。
数式の表示にはMathMLを使っています。また僕は数値計算にあまり詳しくないため、間違いなどありましたら指摘してください。
Durand-Kerner法で多項式の近似解を求めるプログラム
(require (lib "1.ss" "srfi"));SRFI-1のロード
(define PI 3.1415)
(define (square x) (* x x))
(define (dka input) ;Durand-Kerner法
(define formula (if (proper-list? input)
(drop-while zero? input) ;最高次係数が0なら削除
(error "INVALID INPUT")))
(define degree (- (length formula) 1))
(define degree_list (iota degree 1 1))
;次数nに対して、'(1 2 ... n)というリストを作成
(define zcenter
(if (>= 0 degree) (error "INVALID DEGREE")
(- (/ (cadr formula) (* (car formula) degree)))))
(define (get-R list) ;解を囲む円の半径の候補を求める
(let ((denom (car list)))
(map (lambda (x k)
(expt (/ (abs x) denom) (/ 1 k)))
(cdr list) degree_list)))
(define trial-list ;Aberthの初期値
(let* ((angle_base
(/ (* 2 PI (make-rectangular 0 1)) degree))
(R (* degree
(fold max 0 (get-R (kumitate-exp formula zcenter))))))
(map (lambda (x)
(+ zcenter (* R (exp (* angle_base (- x 0.75))))))
degree_list)))
(define (dka-iter list) ;反復
(let ((deriv-list
(map (lambda (x)
(subst (get-form (kumitate formula x)) x))
list))) ;それぞれの近似値における導関数値を求める
(cond ((< (abs-list
(map (lambda (x) (subst formula x)) list))
0.0001) list);それぞれの近似値を与式に代入し、誤差を調べる
(#t (dka-iter
(map (lambda (x y)
(if (= y 0) x (- x (/ (subst formula x) y))))
list deriv-list))))))
(cond ((= degree 1) zcenter)
(#t (dka-iter trial-list))))
(define (kumitate-exp list x) ;組み立て除法を用いて式を変形
(define (kumitate-t-iter list ans)
(cond ((null? list) ans)
(#t (let ((form (kumitate list x)))
(kumitate-t-iter (get-form form) (cons (last form) ans))))))
(kumitate-t-iter list '()))
(define (kumitate list x);組み立て除法
(define (kumitate-iter list v)
(cond ((null? list) '())
(#t (cons (+ (car list) (* v x))
(kumitate-iter (cdr list) (+ (car list) (* v x)) )))))
(cons (car list) (kumitate-iter (cdr list) (car list))))
(define (get-form list) ;組み立て除法後のリストlistから、式の部分を取得
(drop-right list 1))
(define (abs-list list)
(fold + 0 (map (lambda (x) (+ (magnitude x))) list)))
(define (subst form v) ;式formに値vを代入した結果を返す
(last (kumitate form v)))
リスト操作ライブラリ
SRFI-1を使う。proper-list?、 last、drop-while、drop-right、fold、iotaあたりがそう。便利すぎ。
式の表現
多項式は最高次係数を先頭とするリストであらわすことにした。つまりは'(1 -5 6)。
式にたいして値が代入できなきゃどうしようもないので、それを行う関数substを定義。計算には組み立て除法を用いる。
> (subst '(1 -5 6) 2) 0 > (subst '(1 -5 6) 4) 2
組み立て除法
組み立て除法というのはある多項式を (x-a) といった一次式で割るのに適した方法。それを使って商関数Q(x)と剰余Rを求めると、 より、
> (kumitate '(1 -4 -2 34 2 4 3) 2) (1 -2 -6 22 46 96 195) > (kumitate '(1 -2 -6 22 46 96) 2) (1 0 -6 10 66 228)
組み立て除法によって導関数値を求める
組み立て除法を繰り返すと、 ...と更なる変形が得られる。これを使うとaに対するm次の導関数値が比較的簡単に求められる。
いま、f(x)が と変形できたとして、点 x = a における m 次導関数値を求めたい。m 回微分すると次数がm-1以下の項はすべて消えることから、 としてよい。
を考える。これを展開したとき、( x-a )が含まれる項はxにaを代入したとき全て0になるので、結局 が言える。まーここでは一次導関数値しか使わないけど。
Aberthの初期値
ある代数方程式の解が含まれる領域を推定する方法 。直感的には、解の分布の中心zcenterから、その円内にすべての解が含まれるように半径Rの円を描いて、その円周上に初期値を配置する。こうすることで近似値それぞれが真値に近づくスピードをそろえることができる。
解の中心は簡単に求められる。なぜなら、
半径を求めるには、与式を解の中心zcenterで組み立て除法を用いて変形し、r = ( z - zcenter ) と置換することで円の半径 r が解となるような方程式を立てればよい。
結局、n個の初期値を円周上に等間隔に配置するには を j = 1～n について計算すればよい。
こうして求めた初期値はリストとしてtrial-listに格納される。
(define (get-R list) ;解を囲む円の半径の候補を求める
(let ((denom (car list)))
(map (lambda (x k)
(expt (/ (abs x) denom) (/ 1 k)))
(cdr list) degree_list)))
(define trial-list ;Aberthの初期値
(let* ((angle_base
(/ (* 2 PI (make-rectangular 0 1)) degree))
(R (* degree
(fold max 0 (get-R (kumitate-exp formula zcenter))))))
(map (lambda (x)
(+ zcenter (* R (exp (* angle_base (- x 0.75))))))
degree_list)))
Durand-Kerner法
上の方法で求めた初期値それぞれについて、ニュートン法による近似を行う。本来は導関数値をさらに近似する。
とすると、その一次導関数は ここで、 とすると、 の含まれる項は全て消えるので、 このとき、もまた言えるので、
反復処理
(define (dka-iter list) ;反復
(let ((deriv-list
(map (lambda (x)
(subst (get-form (kumitate formula x)) x))
list))) ;それぞれの近似値における導関数値を求める
(cond ((< (abs-list
(map (lambda (x) (subst formula x)) list))
0.0001) list);それぞれの近似値を与式に代入し、誤差を調べる
(#t (dka-iter
(map (lambda (x y)
(if (= y 0) x (- x (/ (subst formula x) y))))
list deriv-list))))))
まず組み立て除法を用いて各近似値に対応する導関数値をそれぞれ計算し、deriv-listとする。それぞれの近似値を与式に代入してみて、全てが0に近ければ十分収束しているとして解を出力、終了。まだダメそうならmapを用いて、近似値全てをニュートン法で再近似。
導関数値リストはもともと別の場所でも使うはずだったので一次変数にしてあるが、その部分を消したためdka-iterの反復呼び出しの際に直接計算してもいい。でもそうするとすごいわかりにくくなると思うのでそのまま。上で示したように、点aにおける一階導関数値を求めるにはまずkumitateを用いて一回変形し、変形後のリストからget-formを用いて式の部分のみを抜き出し、それにaを代入すればよい。
結果
について、
> (dka '(1 -2 1)) (1 1) > (dka '(1 -3 -10 24)) (3.9999998318967864+2.8387151406426836e-008i -3.000000000000001+1.441359404959568e-017i 2.0000000001780975+1.4040173748418967e-010i) > (dka '(1 -5 -21 -22 -16)) (8.0+0.0i -0.4999999999911601+0.8660254037577851i -2.0000000003879568+5.235939426694963e-010i -0.500000298422228-0.8660254066797882i) > (dka '(8 24 -25 3 0 98)) (1.151630579359255+0.817203444821391i -0.7099165021391798+1.038229298587666i -3.883428154439853-9.183549615799121e-041i -0.7099165035189984-1.0382292988192343i 1.151630579277084-0.8172034448086737i)
それぞれ誤差は、
0 2.3890323896364474e-006 7.669790640016116e-006 4.0979890195062324e-007
もっとシンプルにできそうな予感はするんだけど思いつかないなー。
でてた。差分がけっこう多いようだったので、バックアップとってからmtフォルダを書き換え。mt-config.cgi、plugins、dbあたりから必要なものを持ってきたらそのまま使えた。前アップデートしたときは書き換え以外に何か必要だったような気がするんだけど。
集合への30講を読み終えました。整列集合あたりからは僕の想像力の限界を突破してるのでほとんど眼で追っただけだけれども、その難しさはわかった気がする。本書のコラムの中に、2つの異なる実数をとるとは、現実に何を意味しているのだろうか。という言葉があるんですが、うーん、何を意味しているんでしょうね。位相への30講が終わったらもう一回読んでみようと思います。
枝葉の部分でちょっと気になったのは、この本は連続体仮説に連続体「仮設」って文字をあててるんですよね。英語がCH、Continuum Hypothesisだから訳としては仮説なんだろうけど(Googleなんかでも仮説が圧倒的に多い)、証明も反証もできないという意味を考えて仮設にしてあるのかな。
ブラウザで表示するには
まずはwindows。MozzilaやFirefoxの場合、フォントをインストールするだけでよい。MathML at MITに置いてあるMath Fonts to Install for Mozilla/Netscapeを使えば必要なフォントをまとめてインストールすることができる。
入力するには
MathMLの書式を覚えるのも手打ちするのもちょっと面倒なので、入力にはTeX2MathMLを使わせてもらう。指定されたウィンドウにTex形式を入力すると、リアルタイムで数式をレンダリングするとともにMathMLに変換してくれる。Texがそもそもわかんないんですな人はこのへん。
DOCTYPE宣言
こんな感じ。
<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.1 plus MathML 2.0//EN"
"http://www.w3.org/TR/MathML2/dtd/xhtml-math11-f.dtd" [
<!ENTITY mathml "http://www.w3.org/1998/Math/MathML">
]>
<html xmlns="http://www.w3.org/1999/xhtml" xml:lang="ja">
IEだと]>の部分が見えてるなー。でもMozillaではこうしろって書いてあったのでそっち優先。
MIMEの設定
せっかくXHTML1.1で書いたんだし、拡張子をxhtmlにしてMIMEはapplication/xml+xhtmlで送ってやるのがよいと思うのだが、これをやるとIEで表示ができずファイルのダウンロード画面になってしまう(ユーザー側でレジストリをいじれば表示できるようだが)。
text/xmlで送るとW3C側のDTDがおかしいらしく、IEで
'http://www.w3.org/TR/xhtml11/DTD/xhtml11.dtd'の実行エラーです。ライン 85、位置 2
なんてエラーがでる。Firefoxでエラーが出ないのはそもそもDTDのチェックをしてないかららしい……それでいいのか。
text/htmlならIE、Firefoxどちらでも閲覧はできるけれども、MathMLが綺麗に表示されない。どうしたもんかな。
コンテントネゴシエーション
コンテントネゴシエーションについてはこちら。.htaccessを以下のように書き換え。さくらだとユーザー側でOptionsは使えないが、MultiViewsはあらかじめ指定してくれてるみたいなので下二行だけでよい。
Options +MultiViews AddType text/html html AddType application/xhtml+xml xhtml
これでindex.htmlへのリクエストに対し、IEへはindex.html.html、Firefoxへはindex.html.xhtmlを返すことができる。あとはファイルを準備するだけ。
MovableTypeの設定
当たり前だがバックアップとってから。メインページは使ってるテンプレートの出力ファイル名をindex.html.htmlに、テンプレートをもうひとつコピーしてそっちはindex.html.xhtmlにすればOK。
アーカイブは新規でマッピングを作る。アーカイブの拡張子はhtmlのままで、アーカイブ・マッピング側の出力フォーマットの末尾に.htmlと.xhtmlを足したものをそれぞれ作ってやればよい。
こんな感じ。一番上はBlogをFlash化しようとしたときの名残なので気にしないでください。今まで使ってたindex.htmlが残ってるとそっちが表示されてしまうので削除。
リンクの修正
上のようにマッピングを作ると、優先ラジオボタンにチェックのついたURLにMovableType側がリンクを貼ってしまう。これを直したい。
Permalinkの修正
remove-ext-from-permalink.plというPermalinkの拡張子をとっぱらうプラグインがあったので使わせてもらう。これは本来もっとしっかりしたURIに対して使うものみたいで、index.html.*のようなURLじゃないと書き換えてくれない。これをそれ以外のファイル名に対して使うには、ちょっと正規表現を書き換えてやればいい。
ヒント : (.*?)と\1
Archivelinkの修正
これも上と同じようにやろうと思ったんだけど、プラグインの中でどう指定していいのかわからなかったので力技にでる。カテゴリを並べ替えるのに使わせてもらっている、文字列の最初n文字を削るプラグインcutfirstcharを書き換えて、最後のn文字を削除するようにする。上のマッピングの場合.xhtmlへのリンクが張られているので、最後6文字を切り取ってやればよい。
ヒント :
> $string = substr('ABCD', 0, 3);
> print $string;
ABC
> print length($string);
3
そんな感じで
FirefoxならmathMLも綺麗に表示されるし、結果としてxhtml化もできて満足です。今までは個別エントリーの書き出しのことなんてなんも気にせずにMovableTypeのデフォルトのまま行っていたんだけれども、クールなURIについてもちょっと考えさせられたし。つってもこんなことやってないでまずはもっと整ったテンプレートを作れって話だな。
数学を勉強しています。集合・位相に手をつけようと思い、とりあえず図書館へ行って志賀浩二先生の「集合への30講」、「位相への30講」二冊を借りてきた。アレフってどこかでやったな、と思って昔の教科書を調べると、計算機科学入門(アービブ)という本に8Pくらい載ってた(同名でゴールドシュレーガーという人が書いてるのもある。こっちはハードウェアや人工知能にも触れる内容)。どうやら情報基礎学という講義で習ったことがあるようだ。読み返してみるとこの本けっこう盛りだくさんにいろいろ載ってる。まーそんなことはどうでもよくて、集合への30講を読み始めましたよ、という話。
空でない集合Γから、任意の集合Aに対して写像
が与えられたとき、Γを
ある中学校の生徒全員の集合をAとして、学年{1, 2, 3}をΓとすれば、 によってある1学年の生徒の集合が取り出せる。
このとき、その直積集合は と書き、上の例では1学年からひとりずつ生徒を選ぶ組み合わせを意味する。
The Java Programming Language, Fourth Editionをちびちび読んでいる。さっきChapter6 Enumeration Typesを読み終えた。enumの各要素にはメソッドが定義できるので、Stateパターンぽいことができる。
Time.java
public enum Time {
MORNING {
public void greet() {
System.out.println("おはよう");
}
public Time next() {
return DAYTIME;
}
},
DAYTIME {
public void greet() {
System.out.println("こんにちは");
}
public Time next() {
return NIGHT;
}
},
NIGHT {
public void greet() {
System.out.println("こんばんは");
}
public Time next() {
return MORNING;
}
};
public abstract void greet();
public abstract Time next();
}
Test.java
public class Test {
public static void main(String[] args) {
Time t = Time.MORNING;
t.greet();
t = t.next();
t.greet();
t = t.next();
t.greet();
t = t.next();
t.greet();
}
}
コンストラクタも使うとこんな感じになる。greet()は各enum要素内でオーバーライドすることもできるようだ。
Time.java
public enum Time {
MORNING("おはよう") {
public Time next(){
return DAYTIME;
}
},
DAYTIME("こんにちは") {
public Time next(){
return NIGHT;
}
},
NIGHT("こんばんは") {
public Time next(){
return MORNING;
}
};
String message;
private Time(String message) {
this.message = message;
}
public void greet() {
System.out.println(message);
}
public abstract Time next();
}
学校の課題プログラムくらいならこれで十分だろう。Stateパターンを使えってところでこっちを使った場合、どういうことになるかは知らないが。
そういえば先日元バイト先の人と飲む機会があったのだが、情報系の専門学校の子が「javaをフローチャートで習ってる」と言っていた。僕はまず自分の耳とそして彼女の言葉を疑い、それってシーケンス図のことじゃないの、と何度も確認したがどうやら本当にフローチャートらしい。フローチャートって情報処理試験以外でも使われてるんだ、と少し感動した。
