図解・セグメント木
ARC045#B がセグメント木で殴るタイプの問題だった(頭が良い人は 2 回 imos ると簡単に解ける)んですが僕は頭が悪いのでそろそろ自分でセグメント木を実装できるようになろうと思い立ったのでした.が,調べてみたところいろんなサイトの説明を読んでもなかなか理解できず…….複数のページを見比べて数時間かけてようやく実装できたのですが,どうも文章でしか解説されてないせいで全然理解できなかった気がするのでここではちゃんと図(のようなもの)で示しながら説明していきたいと思います.
セグメント木とは
「配列に対する操作の結果をセグメントごとに保持したもの」を木にしたもの.例えば [2, 3, 8, 4, 1, 5, 7, 6]
という配列があったとして,「l 番目から r 番目の要素の中で最小の値は何?」というクエリが何回か投げられるとする (Range Minimum Query: RMQ).これぐらいのサイズの配列なら単純に全部見比べればいいけど,例えば数が 10 億個並んだ配列の中で「1 番目から 10 億番目までの中で……」みたいなことを 10 万回言われると絶対 TLE で落ちる.で,「あらかじめ「全体の中で最小」「右半分で最小」「左半分で最小」……みたいに再帰的に 2 分割しながら結果を持っとけばいいよね」というアイデアなのがセグメント木.上の配列を使うと,
[ 1 ] 0 [ 2 ][ 1 ] 1, 2 [ 2 ][ 4 ][ 1 ][ 6 ] 3, 4, 5, 6 [2, 3, 8, 4, 1, 5, 7, 6] 7, 8, 9, 10, 11, 12, 13, 14 (k) 0, 1, 2, 3, 4, 5, 6, 7 (i) (n = 8)
これは完全 2 分木になっていて,根ノード=min(左の子,右の子) になっていることが分かる.さらに完全 2 分木にすると何が嬉しいかというと,これを 1 次元の配列で扱えるようになる.上図右は左の 2 分木をルートから BFS で 0-indexed の添字を付けたもので,左の子=2 * 親 + 1,右の子=2 * 親 + 2になっている.また,元の配列の長さが 8 で,元の配列の 0 番→ 2 分木配列の 7 番,1 番→ 8番,……という関係なので,元の配列の長さ (8) を n とすると,元の配列の i 番目の要素は 2 分木配列では k = i + n - 1 番目になっている.これで 2 分木配列と元の配列の間でインデックスの変換ができるようになった.
値の更新
セグメント木で出来るのは配列内の 1 要素の更新と,木のある範囲に対する reduce 操作である(max, min, gcd, sum など).要は親=操作(左の子,右の子) が常に成り立てばいい.完全 2 分木では,k 番目のノードの親ノードの番号は (k - 1) / 2 (切り捨て)で求まる.ということで,値の更新の擬似コードは以下のようになる:
def update(元の配列の番号, 新しい値) 2分木内の番号 k = 元の配列の番号 + n - 1 2分木配列[k] = 新しい値 while k > 0: # 1個上の親に上がる k = (k - 1) / 2 # 子を使って親の値を更新する 2分木配列[k] = 操作(2分木配列[2 * k + 1], 2分木配列[2 * k + 2])
値の取得
[ 1 ] 0 [ 2 ][ 1 ] 1, 2 [ 2 ][ 4 ][ 1 ][ 6 ] 3, 4, 5, 6 [2, 3, 8, 4, 1, 5, 7, 6] 7, 8, 9, 10, 11, 12, 13, 14 (k) 0, 1, 2, 3, 4, 5, 6, 7 (i) (n = 8)
例えば,上図で「元の配列で言う 2 番から 7 番までの中で最小の値が欲しい」と言われたら,2 分木配列の添字を見て,直感的に「2 分木配列で言う 2 番,4 番のうち小さい方を返せば良さそう」と思うだろう.思うことにしよう.で,これはどういう直感に基づいているかというと「2 分木配列のノードの中で,欲しい範囲に完全に包含されておりかつもっとも大きい範囲を担当するノード」を返せばいいことになる(これはたぶん iwi さんのスライド 53 枚目が一番分かりやすい).そのためには,2 分木配列の k 番目のノードが元の配列で言う何番から何番までを担当しているかを教えながら再帰的に 2 分探索する必要がある.ということで擬似コードは以下の通り:
# この補助関数を getValue(a, b) = __getValue(a, b, 0, 0, n - 1) と呼べば # 上から順に探索できる def __getValue(元の配列の左端 a, 同右端 b, 2分木配列中の番号 k, k番が担当している元の配列の左端 l, 同右端 r): if k番が担当している元の配列の範囲 l:r が a:b に完全に含まれている: return 2分木配列[k] if l:r が a:b から完全に外れている: return 単位元(minならint.MAX,sumなら0とか) return 操作(__getValue(a, b, 2 * k + 1, l, (l + r) / 2), __getValue(a, b, 2 * k + 2, (l + r) / 2 + 1, r))
この関数では,特に l と r が 2 分木内の添字ではなく元の配列の添字を指していることに注意.また,他の例ではよく半開区間で考えると左の子の右端と右の子の左端が同じになって分かりやすいみたいなことが書いてあったけど,僕にはそれをプログラム内でどう扱えばいいのか分からなかったので普通に閉区間にした.ぶっちゃけ開区間とか閉区間の概念を理解していればどっちでもいいと思う.
プログラム例を挙げるだけでは再帰の処理が何をやっているのか分かりづらいので,例を出して処理を追ってみよう.例えば,getValue(2, 7) とすると補助関数の中では変数が以下のように動いている.
0 1, 2 3, 4, 5, 6 7, 8, 9, 10, 11, 12, 13, 14 (k) 0, 1, 2, 3, 4, 5, 6, 7 (i) (n = 8) 2, 3, 8, 4, 1, 5, 7, 6 (元の配列) a b l r k = 0 => min(4, 1) => 1 l r k = 1 => min(MAX, 4) => 4 l r k = 3 => MAX l r k = 4 => 4 l r k = 2 => 1
2 分木配列の初期化
2 分木配列を初期化するときは,元の配列の長さ以上で最も小さい 2 の冪を求めたあと(完全 2 分木は末端のノード=葉の数が必ず 2 の冪になる),それ * 2 - 1 の長さで中身を単位元で埋めた配列を用意して,k = i + n - 1 番目から順に update していけばいい.たぶんコンストラクタに元の配列を渡してコンストラクタの中で全部やればいいと思う.
まとめ
- セグメント木は,大きい配列の中の部分配列に対する reduce 操作を完全 2 分木を使って出来るようにしたもの
- 値を更新するときは 2 分木を下から上に辿り,値を求めるときは 2 分木を上から分割しながら辿っていけばいい
間違ったこと・分かりづらいことを言ってるところがあったらコメントください.随時修正します.
ちなみに
これを Python 3 で書いて提出したら TLE になった.悲しい.
Submission #525117 - AtCoder Regular Contest 045 | AtCoder
参考
www.slideshare.net
Codeforces Round 313 Div. 2 AからC
Python 3で参加したらデバッグが大変だった.公式のEditorialとほとんど同じだけど一応…….
560A Currency System in Geraldion
貨幣が5種類与えられるので,その貨幣をそれぞれ0回以上使った時にどうしても表現できない最小の値を答える.任意の値を表現できる場合は -1 を出力する.
解法
1 があったら 1 を無限に使うことで任意の値を表現できる.1 が無かったら他の何を使っても 1 を表現できない.これは最小の unfortunate sum.
プログラム
さすがにここまで単純とは思わなくて怖かったけど普通に通った.ちなみに今見たら他の誰かのコードと一致したのか Skipped になってた.そらこんな簡単なの誰かと一致するだろ.だからレート下がったのか.
n = int(input()) a = map(int, input().split()) if 1 in a: print(-1) else: print(1)
560B Gerald is into Art
長方形の絵を2つ買った.長方形のでかいフレームの中に収まるかどうかを判定したい.
解法
(横に並べる・縦に並べる),(1つ目の回転),(2つ目の回転)で 通りある.左上に詰めて置くと考えればいいので,横に並べる時は横の長さが sum で縦の長さが max,縦に並べる時は縦の長さが sum で横の長さが max になることが分かる.8通り全部やる.
プログラム
愚直も愚直.
a1, b1 = map(int, input().split()) a2, b2 = map(int, input().split()) a3, b3 = map(int, input().split()) if (max(a2, a3) <= a1 and b2 + b3 <= b1) or \ (max(a2, b3) <= a1 and b2 + a3 <= b1) or \ (max(b2, a3) <= a1 and a2 + b3 <= b1) or \ (max(b2, b3) <= a1 and a2 + a3 <= b1) or \ (a2 + a3 <= a1 and max(b2, b3) <= b1) or \ (a2 + b3 <= a1 and max(b2, a3) <= b1) or \ (b2 + a3 <= a1 and max(a2, b3) <= b1) or \ (b2 + b3 <= a1 and max(a2, a3) <= b1): print("YES") else: print("NO")
560C Gerald's Hexagon
正三角形から成る六角形の中に詰まっている三角形の数を数える.
想定解法
始めに,正三角形から成る正三角形は,上の列からそれぞれ個数を数えると 1, 3, 5, … というように奇数列になっていることが分かる.奇数列は で表せるので,この sum を取ると になることが分かる.問題で与えられるような六角形は,1辺が の大きい正三角形から上と右下と左下の小さい正三角形を削り取ったものと考えられるので, で答えが求められる.Editorial もうちょっと詳しく書いてくれ.
プログラム
想定解法とは違うやり方で提出していました.だいたい同じですが,平行四辺形を作って右上と左下を削るみたいな計算をしました.上底の長さは で,1列にはその2倍の三角形が詰まっているので,平行四辺形の中にはそれに高さ を掛けた分だけ三角形が詰まっています.後は右上と左下を上と同じように引くだけ.奇数列の総和が という公式に辿り着くまでの計算を渋って補助のメソッドを書いたのがダサい.ちなみにこの解法の式を整理したら とかいう基本対称式の線形結合みたいな式が出てきて特別綺麗でもなく何とも言えない気持ちになりました.
def f(x): ans = 0 for i in range(1, x + 1): ans += 2 * i - 1 return ans a = list(map(int, input().split())) print((a[0] + a[1]) * 2 * (a[1] + a[2]) - f(a[1]) - f(a[4]))
DとEは愚直に書いたら両方TLEで落ちた.お疲れ様でした.C問題はDiv 1のA問題なわけだしこれからもコンスタントにDiv 2で3完していきたいなあと思いました.
PRML 演習9.3
こんなところにもやるだけの問題が.
混合ガウス分布の問題.2値確率変数 は1-of-K符号化されていて,
とする.このとき だから,
総積の部分は典型的な1-of-K符号化の確率の形になっている.これを のときから のときまでのものを全部足し合わせればいいので, についての総和になる.よって
これ以外の問題は全然分からない.
ARC 041 AとB
この前のACM-ICPC国内予選で(主にコーディング面で)全然貢献できなかったという悔しい思い出を機に競プロモチベが若干上がっている昨今。CFの問題をSolver順に解いたりRoundに参加したりしてますが、久しぶりにARCに真面目に参加しました(WORKING'!!一挙放送を見ながらではありましたが)。
A問題 コインの反転
表向きのコイン( 枚)と裏向きのコイン( 枚)を複数枚渡されて、「ちょうど 枚ひっくり返せ」と言われた時にどれだけたくさんのコインを表向きにできるか。
- のとき:裏向きのコインを 枚ひっくり返せばいい→表向きのコインは 枚になる
- のとき:とりあえず裏向きの 枚は全部ひっくり返す・表向きのコインのうち被害を受けるのは 枚→表向きのコインは 枚
公式解説スライドは絶対値を使って に依らない想定解を出してたけどたぶん場合分けして展開したらこうなるんだと思います(適当)。実際2つ目式の括弧の中の符号を反転すれば1つ目の式になるっぽい。雑魚だからコンテスト中に絶対値の中の正負がとか考えてる余裕はない。
import std.stdio; import std.string; import std.algorithm; import std.conv; import std.range; import std.math; void main() { int x, y; readf("%d %d\n", &x, &y); int k = readln.chomp.to!int; if (k <= y) { writeln(x + k); } else { writeln(x + 2 * y - k); } }
B問題
「わっ!」ってやると四方に分裂するアメーバがいる。「わっ!」って言った後の盤面が与えられるので、「わっ!」って言う前の盤面を出力する。
「あるセルの上下左右全てにアメーバがいたらそのセルにはアメーバがいた」と考えられる。左上から右下に向かって順番に「上下左右のアメーバのうち一番少ない数( とする)をそのセルにいたアメーバの数にし、上下左右から を引く」操作を続けると矛盾なく解ける。
import std.stdio; import std.string; import std.algorithm; import std.conv; import std.range; import std.math; void main() { int n, m; readf("%d %d\n", &n, &m); int[][] b; for (int i = 0; i < n; i++) { b ~= readln.chomp.split("").map!(to!int).array; } int[][] a = new int[][](n, m); for (int i = 1; i < n - 1; i++) { for (int j = 1; j < m - 1; j++) { int up = b[i][j - 1]; int down = b[i][j + 1]; int right = b[i + 1][j]; int left = b[i - 1][j]; if (up > 0 && down > 0 && right > 0 && left > 0) { int min_am = min(up, down, right, left); a[i][j] = min_am; b[i][j - 1] -= min_am; b[i][j + 1] -= min_am; b[i + 1][j] -= min_am; b[i - 1][j] -= min_am; } } } foreach (e; a) { e.map!(to!string).join("").writeln; } }
C
RRR....LLLL的なセグメントに分けるところまでは思い付いたけど「常に多い方を少ない方に寄せればいい」という簡潔な最適解まで思いつけなかった。
D
読んでない。
Sublime Text 3のD言語開発パッケージDKitを使う
概要
Sublime Text 3 で D言語を開発するためのツール DKit を Mac で導入します.結構前に出たパッケージだし紹介記事ぐらいあるだろうと思っていたのですが,UNIX 系の物は無かったので一応.
手順
DKit は D Completion Daemon(通称DCD)と連携して動作します.Homebrew にいるので brew install dcd
で入ります.
README にもありますが,DKit はまだ開発途中のため,Package Control からインストールすることはできません.そのため,Packages ディレクトリに直接展開します.
# Packages ディレクトリの在処はメニューの Sublime Text → Preferences → Browse Packages... から分かります $ cd ~/Library/Application\ Support/Sublime\ Text\ 3/Packages $ git clone https://github.com/yazd/DKit
次に,設定ファイルを変更します.例によって Sublime Text → Preferences → Package Settings → DKit → Setting − User をいじります.Default に設定されている値を上書きするだけですが,DMD を Homebrew でインストールした場合は以下のようになります.
{ "dcd_path": "/usr/local/bin", "include_paths": [ "/usr/local/include/d2/" ] }
これで準備完了です.D 言語のファイルを開いてみます.
おお…….続けて見ていきましょう.
ちゃんと名前空間を考慮して補完されていきます.まあここらへんは全部 DCD から来た値を受け取っているだけだと思いますが.
スニペットも用意されています.
関数もちゃんとインポートされたものの中からのみ補完されます.
ファジィ検索もお手の物.
実行する時は Super + B,D - run を選択すると rdmd で実行するので標準出力を確認できます.
キーバインドで Shift + Super + G が Default の方で用意されており,関数の定義元に一瞬でジャンプできます.当然ローカルで定義した関数にもジャンプできます.
たまに Tab を押しても何も起こらないことがあるのですが,補完したかった場所でもう一度 Tab を押すとちゃんと補完されます.どうしてこういうことが起こるのかよく分かりませんが.また,現状 DCD の方で UFCS に対応していないので,Ruby のように気持ちよくメソッドチェーンするみたいなことはまだ出来ませんが,Sublime Text 2 が出た頃に D 言語のスキーマが用意されていなかったことを思えばずいぶん遠くに来たなあと思います.大学に入ってから Emacs 一辺倒だったのですが,移行しない理由の半分が D 言語が書けないことだったので,ここまで出来るとなるとうっかり移行してしまいそうです…….
AOJ CGL_1_A 射影
大学の図書館に入っていたのでこの本をざっと読んでいました.AOJにある幾何の問題を全然解いていないことに気付いたので着手しました.簡単な問題から.
問題
平面上にある点からある直線に下ろした垂線の足を求める.
解釈
幾何的にも解けるけど線形代数を使って解いてみる.下図の が分かれば より求められる.
は点 を通る直線上にあるため,
と表せる.そのため
である( の正負は の値によって変わる).また, と の内積を考えると
だが,図中の直角三角形に着目すると
であることが分かる.したがって
より である.式(1), (2)より,
よって
だから,
の係数がスカラーであることに気付かないと一見「あれ?」ってなる.とにかくここまで来れば最後の式を実装するだけ.
import std.stdio; import std.string; import std.algorithm; import std.conv; import std.range; import std.numeric; void main() { auto point = readln.chomp.split.map!(to!double).array; auto p1 = point[0..2]; auto p2 = point[2..4]; int q = readln.chomp.to!int; for (int i = 0; i < q; i++) { auto p = readln.chomp.split.map!(to!double).array; double[2] p1p = p[] - p1[]; double[2] p1p2 = p2[] - p1[]; double[2] ph = p1p2[] * dotProduct(p1p, p1p2) / dotProduct(p1p2, p1p2) + p1[]; writefln("%.10f %.10f", ph[0], ph[1]); } }
を求めるには,各点同士の座標を引き算するだけ.D言語には配列のスライスを使ったベクトル演算が可能なため,p[] - p1[]
のようなベクトルの定義に沿った書き方が出来る(p[] += 3
などのようなベクトルに対するスカラー演算も可能です).dotProduct
はstd.numeric
に定義されている.D言語は神.
参考
AtCoder Typical Contest 001 をD言語で解いた(A, B問題)
競プロ界隈重要アルゴリズムが大会形式でchokudaiさんのスライドを見ながら提出できるという神コンテストが開催されたので,勉強不足な僕も参加してみました.RubyでDFSを書き始めたらデバッグが全く進まなくて無事死亡したので後日D言語で改めて解いたものをブログに残しておきます.
A問題
迷路をスタートからゴールまで行く.深さ優先探索という問題なので一応それで.迷路と言われたら普通最短解で幅優先な気がするけどTwitter見てたらUnion Findでも解けるみたいな意見を見て「??????」ってなってた.やり方が全く思いつかない.
import std.stdio; import std.range; import std.string; void main() { int H, W; readf("%d %d\n", &H, &W); string[] c; int[] s = new int[2]; foreach (int i, string line; stdin.lines) { c ~= line.chomp; foreach (int j, ch; line) { if (ch == 's') s = [ i, j ]; } } /** * ある地点に訪問済みかどうかを記録しておく配列 * 記録しておかないと同じ地点を永遠にさまよい続ける */ bool[][] visit = new bool[][](H, W); if (dfs(H, W, s[0], s[1], c, visit)) writeln("Yes"); else writeln("No"); } /** * 深さ優先探索で迷路を探索する * H: 高さ, W: 幅, h: 今いる場所の縦位置, w: 今いる場所の横位置 * c: 迷路 v: 記録シート */ bool dfs(int H, int W, int h, int w, string[] c, bool[][] v) { // 迷路からはみ出たらやめる // これを一番最初に弾いとかないと Range violation で実行時に落ちる if (h < 0 || H <= h || w < 0 || W <= w) return false; // 訪問済みだったらやめる if (v[h][w]) return false; // 壁だったらやめる if (c[h][w] == '#') return false; // ゴールに辿り着いたら成功 if (c[h][w] == 'g') return true; // 今たどり着いたマスをチェックしておく v[h][w] = true; // 上下左右のうちどれかがゴール = trueだったらいいので論理和を返す return dfs(H, W, h - 1, w, c, v) || dfs(H, W, h + 1, w, c, v) || dfs(H, W, h, w - 1, c, v) || dfs(H, W, h, w + 1, c, v); }
readf
で改行を読み込み忘れたり,配列のサイズをコンパイル時に確保しといたのに実行時に ~=
で追加して前の方に空白の要素が出来ちゃったりいろいろな凡ミスをしてしまった.いい加減流れを覚えたい.あと最初はスタックを使って書いてみたけど配列を適当に弄ったのでたぶんスタック配列からの取り出しとか追加に時間がかかってTLEで落ちた.関数の再帰呼出しに書き直したら普通に通った.でもスタックオーバーフローで落ちるのこわい.
B問題
N個ある要素をつなげたり,与えられた2つの要素が繋がっているかどうかの判定をしたりしたい.Union Findって今まで名前は知ってたけどどういうものか分からなかったので,今回のお陰で理論と実装が分かってよかった.応用の仕方は全く分からないので良い問題を教えてほしい.
import std.stdio; import std.range; void main() { int N, Q; readf("%d %d\n", &N, &Q); int[] uf = new int[](N); // 最初は各要素が自分を見るようにする // ref を付けないとコピーしたのを見てしまうので意味が無い foreach (int i, ref e; uf) { e = i; } foreach (q; iota(0, Q, 1)) { int p, a, b; readf("%d %d %d\n", &p, &a, &b); if (p == 0) { unite(a, b, uf); } else { if (find(a, b, uf)) writeln("Yes"); else writeln("No"); } } } bool find(int a, int b, int[] uf) { return (root(a, uf) == root(b, uf)); } // 引数で受け取ったUF木を直接いじりたいのでrefを付ける void unite(int a, int b, ref int[] uf) { int x = root(a, uf); int y = root(b, uf); // もし根が違ったら根どうしを繋げる if (x != y) uf[x] = y; } int root(int x, int[] uf) { // 自分自身を参照しているところが根 if (uf[x] == x) return x; else // 代入しながら再帰することで,途中経過で参照したノードを全て // 根に繋げることができる // 例:uf[2] = uf[5] = uf[3] = uf[4] = 4 // で,2と5と3がこれ以降直接4を見てくれるようになる return uf[x] = root(uf[x], uf); }
C問題は脳が途中で理解しようとする努力をやめてしまいました.南無三.