Pythonの練習
Pythonの練習として、PriorityQueueとUnionFindを実装してみた。まだPythonをよくわかってないので間違ったことをしてるかも。
あとテストにunittest使ってみた。
#!/usr/bin/env python # -*- coding: utf-8 -*- class PriorityQueue: """ 優先度付キュー 演算子"<"で比較できるオブジェクトなら何でも格納できる 小さいオブジェクトから順にpopされる """ def __init__(self): """ 空のキューを生成する O(1) """ self.heap = [] def push(self, x): """ xをキューに追加する O(log N) """ self.heap.append(x) self._siftUp() def pop(self): """ キューから一番小さいオブジェクトを取り除いて返す O(log N) """ assert not self.empty() if len(self.heap) == 1: return self.heap.pop() else: x = self.heap[0] self.heap[0] = self.heap.pop() self._siftDown() return x def empty(self): """ キューが空なら真を返す O(1) """ return len(self.heap) == 0 def _siftUp(self): heap = self.heap i = len(heap) - 1 while i > 0: parent = (i - 1) / 2 if parent >= 0 and heap[i] < heap[parent]: heap[i], heap[parent] = heap[parent], heap[i] i = parent else: break def _siftDown(self): heap = self.heap i = 0 while i < len(heap): child1 = i * 2 + 1 child2 = i * 2 + 2 if child1 < len(heap) and child2 < len(heap): if not (heap[child1] < heap[i]) and not (heap[child2] < heap[i]): break elif heap[child1] < heap[child2]: heap[i], heap[child1] = heap[child1], heap[i] i = child1 else: heap[i], heap[child2] = heap[child2], heap[i] i = child2 elif child1 < len(heap) and heap[child1] < heap[i]: heap[i], heap[child1] = heap[child1], heap[i] i = child1 else: break if __name__ == '__main__': import unittest class PriorityQueueTest(unittest.TestCase): def setUp(self): self.q = PriorityQueue() def test(self): q = self.q assert q.empty() q.push(10) assert not q.empty() assert q.pop() == 10 assert q.empty() q.push(10) q.push(5) q.push(100) q.push(20) q.push(77) q.push(0) q.push(140) q.push(123) assert q.pop() == 0 assert q.pop() == 5 assert q.pop() == 10 assert q.pop() == 20 assert q.pop() == 77 q.push(8) q.push(130) q.push(130) q.push(200) assert q.pop() == 8 assert q.pop() == 100 assert q.pop() == 123 assert q.pop() == 130 assert q.pop() == 130 assert q.pop() == 140 q.push(-1000) assert q.pop() == -1000 assert q.pop() == 200 assert q.empty() suite = unittest.TestLoader().loadTestsFromTestCase(PriorityQueueTest) _unit = unittest.TextTestRunner(verbosity=2).run(suite) print _unit
#!/usr/bin/env python # -*- coding: utf-8 -*- class DisjointSets: """ 互いに素な集合族に対して，次の2つの操作を高速に行うデータ構造 - union: 2つの集合を合併する - find: 2つの要素が同じ集合に属しているか調べる """ def __init__(self, n): self.parent = n * [-1] def union(self, x, y): x = self.root(x) y = self.root(y) if x != y: if self.parent[x] < self.parent[y]: self.parent[x] += self.parent[y] self.parent[y] = x else: self.parent[y] += self.parent[x] self.parent[x] = y return x != y def find(self, x, y): return self.root(x) == self.root(y) def root(self, x): y = x while self.parent[x] >= 0: x = self.parent[x] while self.parent[y] >= 0: self.parent[y], y = x, self.parent[y] return x if __name__ == '__main__': import unittest class DisjointSetsTest(unittest.TestCase): def test(self): ds = DisjointSets(5) assert not ds.find(0, 1) assert not ds.find(1, 2) assert not ds.find(0, 3) assert not ds.find(3, 4) ds.union(1, 2) assert ds.find(1, 2) assert not ds.find(0, 1) assert not ds.find(3, 4) ds.union(4, 2) assert ds.find(4, 2) assert ds.find(1, 4) assert not ds.find(0, 1) assert not ds.find(3, 4) ds.union(0, 3) assert not ds.find(0, 1) assert not ds.find(3, 4) assert ds.find(0, 3) ds.union(3, 1) assert ds.find(3, 1) assert ds.find(4, 0) suite = unittest.TestLoader().loadTestsFromTestCase(DisjointSetsTest) _unit = unittest.TextTestRunner(verbosity=2).run(suite) print _unit
Google Code Jam Qualification Round 2010
Google Code Jam Qualification Round 2010に参加した。3問とも通ったらしい。
A. Snapper Chain
結局、T-FFをN個連結して、クロックをK回入れたらどうなりますかっていう話なので、英語さえ読めればコードは楽勝。
#include <cstdio> using namespace std; int main() { int T; scanf("%d", &T); for (int t = 0; t < T; ++t) { printf("Case #%d: ", t+1); int N, K; scanf("%d%d", &N, &K); if (K % (1 << N) == (1 << N) - 1) puts("ON"); else puts("OFF"); } }
B. Fair Warning
一番悩んだ。
条件式を変形していくと
ti + y = 0 (mod T)
ti = -y (mod T)
t1 = t2 = ... = tn = -y (mod T)
t1 = t2 かつ t2 = t3 かつ ... かつ tn-1 = tn (mod T)
t1 - t2 = 0 かつ t2 - t3 = 0 かつ ... かつ tn-1 - tn = 0 (mod T)
となる。最後の式は、Tが{ti - ti+1}の公約数であることを意味する。
逆に、Tが{ti - ti+1}の公約数であれば、条件を満たすyが存在する。(例えば y = (T - (t1 % T)) % T があり、実はこれが最小のy)
よってTの最大値は{ti - ti+1}の最大公約数。
require 'mathn' gets.to_i.times do |c| n, *ts = gets.split.map{|s| s.to_i} g = (ts[1] - ts[0]).abs for i in (1..n-2) g = g.gcd((ts[i+1] - ts[i]).abs) end puts "Case ##{c + 1}: #{(g - ts[0] % g) % g}" end
C. Theme Park
周期性を利用する。
もし、ある2つの時点で、コースターに乗る先頭のグループが同じなら、それ以降は全く同じ状況が起こる。先頭の選び方はN通りしかないので、高々N+1回でコースターに乗る先頭のグループが同じという状況が起こる。その後は、ループが１週したときの儲け * ループする回数 で一気に計算できる(余りの部分は適当になんとかする)。
#include <cstdio> #include <cstring> #include <cassert> using namespace std; long long g[1000]; int turn[1000]; long long money[1000]; int main() { int T; scanf("%d", &T); for (int t = 0; t < T; ++t) { memset(turn, -1, sizeof(turn)); memset(money, -1, sizeof(money)); int R, k, N; scanf("%d%d%d", &R, &k, &N); for (int i = 0; i < N; ++i) { scanf("%lld", g+i); } int head = 0; long long m = 0; for (int i = 0; i < R; ) { if (turn[head] != -1) { long long period = i - turn[head]; long long mouke = m - money[head]; long long times = (R - i) / period; assert(mouke * times >= 0); assert(period * times >= 0); m += mouke * times; i += period * times; memset(turn, -1, sizeof(turn)); memset(money, -1, sizeof(money)); } else { money[head] = m; turn[head] = i; long long n_ride = 0, j; for (j = head; ; j = (j + 1) % N) { if (j == head && n_ride > 0) break; if (n_ride + g[j] > k) break; n_ride += g[j]; } m += n_ride; head = j; ++i; } assert(m >= 0); } printf("Case #%d: %lld\n", t+1, m); } }
SRM 463
Div2 Easy: BunnyPuzzle
うさぎが条件を満たしつつ移動できる場合の数を求める問題。
class BunnyPuzzle { public: int theCount(vector <int> bunnies) { int n = bunnies.size(); int result = 0; for (int i = 0; i < n; ++i) { for (int j = i-1; j >= 0; --j) { int l = bunnies[i] - bunnies[j]; int count = 0; for (int k = i-1; k >= 0 && bunnies[k] >= bunnies[i]-2*l; --k) ++count; if (count == 1) ++result; } for (int j = i+1; j < n; ++j) { int l = bunnies[j] - bunnies[i]; int count = 0; for (int k = i+1; k < n && bunnies[k] <= bunnies[i]+2*l; ++k) ++count; if (count == 1) ++result; } } return result; } };
Div2 Medium: RabbitNumbering
n要素の整数列{a_i}が与えられるので, n要素の順列pの中で, すべてのi≦nに対して p_i ≦ a_i であるようなものの数を求めよという問題。
class RabbitNumbering { public: int theCount(vector <int> maxNumber) { int n = maxNumber.size(); sort(maxNumber.begin(), maxNumber.end()); long long result = 1; for (int i = 0; i < n; ++i) { result *= maxNumber[i] - i; result %= 1000000007; } return result; } };
Div2 Hard: Nisoku
n枚の数字の書かれたカードの山が与えられる。カードの中から任意に2枚取り出し, 2枚のカードに書かれている数字の和か積を書き込んだ新しいカードを山に加える。取り出した2枚のカードは捨てる。この作業を繰り返し, 山にあるカードが1枚になったとき, そのカードに書かれている数字の最大値をもとめよ, という問題。
class Nisoku { public: double theMax(vector <double> cards) { int n = cards.size(); double result = 0; sort(cards.begin(), cards.end()); for (int i = 0; i <= n; i += 2) { double v = 1; for (int j = 0; j < i/2; ++j) v *= cards[j] + cards[i-j-1]; for (int j = i; j < n; ++j) v *= cards[j]; result = max(result, v); } return result; } };
この問題のアルゴリズムそのものはそんなに難しくない。
- 偶数iに対して, カードの山を2つの集合 A = {i番目に小さいカードより小さいカード}, B = {Aに入ってないカード} に分ける.
- Aに含まれるカードから順に一番大きいカードと一番小さいカードを取り出し和をとったものの積を計算する。
- Bの含まれるカードの積を計算する。
- AとBで計算した数字の積がその分割でできる最大値なのですべてのiについてこの値を計算すればいい。
アルゴリズムはこれだけだが, これが正しく動くことの証明が結構大変だった。
- てきとうな証明(pdf) http://kmc.jp/~nojima/files/SRM/nisoku.pdf
[追記] 証明に山の分割について書いてなかったので追記。単純にある数字を閾値にして2つに分割するだけで最大値が見つかるのは, 以下の式が成り立つからです。
a≧0 かつ b≦c ならば (a+c)b ≦ (a+b)c
(証明) (a+b)c - (a+c)b = a(c-b) ≧ 0
最近の.vimrc
"===================================================================== " Vim 設定ファイル "===================================================================== "--------------------------------------------------------------------- " 編集 "--------------------------------------------------------------------- set tabstop=8 set shiftwidth=2 set expandtab set smarttab set cindent set backspace=2 set wildmenu set formatoptions+=mMr "--------------------------------------------------------------------- " 画面表示 "--------------------------------------------------------------------- syntax on filetype on filetype indent on filetype plugin on colorscheme inkpot set number set laststatus=2 set cmdheight=1 set showcmd set notitle set scrolloff=5 set nocp incsearch set statusline=%<%f\ %m%=%l/%L,%v\ %r%{'['.(&fenc!=''?&fenc:&enc).']['.&ff.']'} set foldmethod=marker "--------------------------------------------------------------------- " 検索 "--------------------------------------------------------------------- set ignorecase set smartcase "--------------------------------------------------------------------- " ファイル操作 "--------------------------------------------------------------------- set backup set backupdir=$HOME/.vim/backup set autochdir "--------------------------------------------------------------------- " バッファ "--------------------------------------------------------------------- set hidden "--------------------------------------------------------------------- " キーマップ "--------------------------------------------------------------------- " バッファ・タブを切り替える nnoremap <Space>j :bn<CR> nnoremap <Space>k :bp<CR> nnoremap <Space>h :tabn<CR> nnoremap <Space>l :tabp<CR> " 表示行単位で移動できるようにする nnoremap j gj nnoremap k gk vnoremap j gj vnoremap k gk nnoremap <Down> g<Down> nnoremap <Up> g<Up> inoremap <Down> <C-o>gj inoremap <Up> <C-o>gk vnoremap <Down> g<Down> vnoremap <Up> g<Up> "--------------------------------------------------------------------- " 文字コード・改行コードの設定 "--------------------------------------------------------------------- set fileencodings=ucs-bom,utf-8,iso-2022-jp,euc-jp,cp932,utf-16,utf-16le set fileformats=unix,dos,mac "--------------------------------------------------------------------- " その他 "--------------------------------------------------------------------- " 曖昧な幅の文字を全角として認識する set ambiwidth=double " beep音を消す set vb t_vb= " C settings set cinoptions=:0,p0,t0 set cinwords=if,else,while,do,for,switch,case
パターン認識と機械学習
なんか後輩がPythonで実際にコードを書いてるらしいので、自分も最小二乗法のコードを書いてみることにした。
最初はRubyで書く予定だったけど、せっかくだからRを使ってみようとか思い、ググりながらRで書いてみた。
# サインカーブを描く x <- seq(0, 1, 0.01) plot(x, sin(2*pi*x), type="l", xlim=c(0, 1), ylim=c(-1.2, 1.2), ann=F, col=2) # サンプルをとる N <- 10 # 点の数 x <- (1:N)/N y <- sin(2*pi*x) + 0.2*rnorm(N) # サンプルを描画する par(new=T) plot(x, y, xlim=c(0, 1), ylim=c(-1.2, 1.2), ann=F, col=4) M <- 3 # フィットする多項式の次元 # 最小二乗法で多項式の係数を求める A <- matrix(nrow=M+1, ncol=M+1) for (i in 0:M) { for (j in 0:M) { A[i+1,j+1] <- sum(x^(i+j)) } } Y <- matrix(nrow=M+1, ncol=1) for (i in 0:M) { Y[i+1] <- sum(x^i * y) } w = solve(A, Y) # できた多項式を描画する x <- seq(0, 1, 0.005) y <- c() for (i in 1:length(x)) { y[i] <- sum(w * x[i]^(0:M)) } par(new=T) plot(x, y, type="l", xlim=c(0, 1), ylim=c(-1.2, 1.2), ann=F, col=3)
Rは初めてなので、これでいいのかすごく不安だけど、とりあえず以下の出力が得られたので満足。
SRM 2回目
o x x 173.96pt でレーティングが1116 -> 1041となった。相変わらずひどい。
Div2 Easy: Archery
アーチェリーの得点の期待値を計算する問題。
class Archery { public: double expectedPoints(int N, vector <int> ringPoints) { double prob = 0; for (int i = 0; i <= N; ++i) { prob += ringPoints[i] * ((i+1)*(i+1) - i*i); } return prob / ((N+1)*(N+1)); } };
(includeしてる行がやたらいっぱいあってうっとうしいのでincludeの行は省いた)
Div2 Medium: AgeEncoding
0と1だけからなる文字列candlesLineが与えられる。この文字列をb進数と見たときcandlesLineがageと等しくなるようなbを求めよ。もし、そのようなbが存在しないのなら-1を、複数存在するのなら-2を出力せよ。という問題。
age=1, candlesLine="111111111111111111111111" のような場合は-1を出力するのが正解だけど0を出力してて通らなかった。
以下、本番後に書いたコード。
class AgeEncoding { public: double calc(string candlesLine, double b) { double value = 0; for (int i = 0; i < candlesLine.size(); ++i) { value += (candlesLine[i] - '0') * pow(b, (double)(candlesLine.size() - i - 1)); } return value; } double getRadix(int age, string candlesLine) { for (int i = 0; i < candlesLine.size()-1; ++i) if (candlesLine[i] == '1') goto ok; if (candlesLine[candlesLine.size()-1] == '1' && age == 1) return -2; else return -1; ok: double lo = 0, hi = 101; for (int k = 0; fabs(hi - lo) > 1e-10 && k < 1000; ++k) { double mid = (lo + hi) / 2; double value = calc(candlesLine, mid); if (value < age) { lo = mid; } else { hi = mid; } } return lo <= 100 && lo >= 1e-9 ? lo : -1; } };
Div2 Hard: SteeplechaseTrack
有向グラフが与えられるので、頂点数がN以下のウォークのうち最大の長さをもつウォークの長さを答えよ。という問題。(正確にはちょっと違うけど)
Bellman Ford法のような動的計画法で解こうとしてたけど時間が足りなかった。
#define REP(i, n) for (int i = 0; i < (int)(n); ++i) #define FOR(i, a, b) for (int i = (int)(a); i < (int)(b); ++i) const int INF = 100000; class SteeplechaseTrack { public: int maxComplexity(vector <string> fences, vector <string> tracks, int N) { REP(i, tracks.size()) REP(j, tracks[i].size()) tracks[i][j] -= '0'; REP(i, fences.size()) REP(j, fences[i].size()) fences[i][j] -= '0'; int n = fences.size(); vector<vector<int> > dp(N, vector<int>(n, -INF)); REP(i, n) if (fences[i][1] != 0) dp[0][i] = fences[i][0] + fences[i][1]; FOR(k, 1, N) REP(i, n) REP(j, n) if (tracks[i][j] != 0) dp[k][j] = max(dp[k][j], dp[k-1][i] + fences[j][0] + tracks[i][j]); int result = -INF; REP(k, N) REP(i, n) if (fences[i][2] != 0) result = max(result, dp[k][i] + fences[i][2]); return result < 0 ? -1 : result; } };
初めてSRMに参加した
初めてSRMに参加しました。結果は o x x で 230.12pt 510位 という残念なものでした。次回に向けて今回の問題について書いておきます。(ソースコードは本番後に清書したもの)
Div2 Easy: TrappingRabbit
一番近い罠までのマンハッタン距離を求める問題。普通に求めて終了。この問題はできた。
#include <vector> #include <algorithm> using namespace std; class TrappingRabbit { public: int findMinimumTime(vector<int> trapX, vector<int> trapY) { int result = 100000; for (int i = 0; i < trapX.size(); ++i) result = min(result, abs(trapX[i] - 1) + abs(trapY[i] - 1)); return result; } };
Div2 Medium: ColoringRectangle
長方形が一つと、赤と青の円盤が何枚か与えられるので、ある条件を満たしながら、長方形全体を円盤で覆うのに必要な円盤の最低枚数を求めよという問題。
条件から円盤は半径が大きい順に使えばいいことが分かり、さらに条件から赤と青を交互に使わないといけないので、結局、最初に赤い円盤を使う場合と青い円盤を使う場合の2通りしかないことが分かる。よって両方の場合について計算してminをとればいい。
本番では条件を一部誤解していて撃墜されてしまった。
#include <cmath> #include <vector> #include <algorithm> #include <functional> using namespace std; const double eps = 1e-5; class ColoringRectangle { public: int attempt(int width, int height, vector<int>& disk1, vector<int>& disk2) { vector<int> disk[2] = { disk1, disk2 }; double dist = 0; for (int i = 0; ; ++i) { int j = i >> 1; int k = i & 1; if (j >= disk[k].size()) break; dist += sqrt(disk[k][j]*disk[k][j] - height*height); if (dist >= width - eps) return i+1; } return -1; } int chooseDisks(int width, int height, vector<int> red, vector<int> blue) { sort(red.begin(), red.end(), greater<int>()); sort(blue.begin(), blue.end(), greater<int>()); return (int)min((unsigned)attempt(width, height, red, blue), (unsigned)attempt(width, height, blue, red)); } };
disk[k][j] < height の場合の処理を書いてないけど、この場合は sqrt(disk[k][j]*disk[k][j] - height*height) がNaNになり、NaNと何かを比較すると必ずfalseになるので、特に処理を書かなくても上手くいく。
Div2 Hard: NameInput
予測文字列と名前が与えられるので、予測文字列を使って名前を入力するのに何ステップいるかという問題。
n = name.size(), m = predictionSequence.size() とすると、インプットをn*mのサイズの地図とみなして幅優先探索すれば、O(nm)で最短経路が求まる。
本番では問題の規模を誤解してて(またかい)、O(nm^2)のDPで解いててTLE.
#include <string> #include <vector> #include <algorithm> #include <queue> #include <cstring> using namespace std; class NameInput { public: int dist[2501][2500]; int countUpDownKeyPresses(vector<string> predictionSequence, vector<string> names) { string pred; for (int i = 0; i < predictionSequence.size(); ++i) pred += predictionSequence[i]; string name; for (int i = 0; i < names.size(); ++i) name += names[i]; memset(dist, -1, sizeof(dist)); queue<pair<int, int> > Q; Q.push(make_pair(0, 0)); dist[0][0] = 0; while (!Q.empty()) { int n = Q.front().first; int p = Q.front().second; Q.pop(); if (n < name.size() && name[n] == pred[p] && dist[n+1][p] == -1) { dist[n+1][p] = dist[n][p]; Q.push(make_pair(n+1, p)); } int up = (p - 1 + pred.size()) % pred.size(); if (dist[n][up] == -1) { dist[n][up] = dist[n][p] + 1; Q.push(make_pair(n, up)); } int down = (p + 1) % pred.size(); if (dist[n][down] == -1) { dist[n][down] = dist[n][p] + 1; Q.push(make_pair(n, down)); } } unsigned result = (unsigned)-1; for (int i = 0; i < pred.size(); ++i) result = min(result, (unsigned)dist[name.size()][i]); return (int)result; } };
ついでにDiv1にも挑戦してみた。Div1のEasyはDiv2のMediumと同じ問題っぽいので省略。Hardは解き方分からなかった。
Div1 Medium: BuildingCities
二次元平面上の街の集合が与えられる。街と街はその間のユークリッド距離がmaxDirect以下の時にその二つの街をつなぐ道が存在する。ここで、街0から街n-1までの最短経路の長さをmaxTravel以下にするように新たな街をいくつか作りたい。最低いくつの街を作ればよいか求めよという問題。
新たに作る街の数kに対してk回のDijkstra法を行うことで解いた。
#include <vector> #include <algorithm> #include <cmath> #include <queue> using namespace std; class BuildingCities { public: double dp[3001][50]; bool visited[3001][50]; int euclid2(int x1, int y1, int x2, int y2) { int dx = x2 - x1; int dy = y2 - y1; return dx*dx + dy*dy; } double euclid(int x1, int y1, int x2, int y2) { return sqrt(euclid2(x1, y1, x2, y2)); } int findMinimumCities(int maxDirect, int maxTravel, vector<int> cityX, vector<int> cityY) { int n = cityX.size(); if (euclid2(cityX[0], cityY[0], cityX[n-1], cityY[n-1]) > maxTravel*maxTravel) return -1; for (int i = 0; i <= 3000; ++i) { fill(visited[i], visited[i] + n, false); fill(dp[i], dp[i] + n, 1e10); } priority_queue<pair<double, int> > Q[3001]; Q[0].push(make_pair(0, 0)); dp[0][0] = 0; for (int k = 0; ; ++k) { while (!Q[k].empty()) { double cost = -Q[k].top().first; int city = Q[k].top().second; Q[k].pop(); if (visited[k][city]) continue; visited[k][city] = true; for (int i = 0; i < n; ++i) { if (!visited[k][i]) { double dist = euclid(cityX[city], cityY[city], cityX[i], cityY[i]); int need = dist <= maxDirect ? 0 : (int)ceil(dist / maxDirect) - 1; double new_cost = cost + dist; if (k+need <= 3000 && new_cost < dp[k+need][i]) { dp[k+need][i] = new_cost; Q[k+need].push(make_pair(-new_cost, i)); } } } } if (dp[k][n-1] <= maxTravel) return k; } } };
