LCAを求める方法はいくつかありますが、今回はダブリングを使う方法とHL分解（Heavy Light Decomposition）を使う方法を紹介します。【追記】オイラーツアーを使う方法を追加しました。【追記】Link-Cut Treeを使う方法を追加しました。
ダブリングによる実装
template <typename T = int> struct LowestCommonAncestor { struct Edge { int to; T cost; Edge(int to, T cost) : to(to), cost(cost) {} }; int n; int l; vector<vector<Edge>> es; vector<vector<int>> par; vector<int> dep; vector<T> dists; LowestCommonAncestor(int n) : n(n), es(n), dep(n), dists(n + 1) { l = 0; while ((1 << l) < n) l++; par.assign(n + 1, vector<int>(l, n)); } void add_edge(int v, int u, T c = 0) { es[v].emplace_back(u, c); es[u].emplace_back(v, c); } void dfs(int v = 0, int d = 0, T c = 0, int p = -1) { if (p != -1) par[v][0] = p; dep[v] = d; dists[v] = c; for (auto e : es[v]) { if (e.to == p) continue; dfs(e.to, d + 1, c + e.cost, v); } } void build() { dfs(); for (int i = 0; i < l - 1; i++) { for (int v = 0; v < n; v++) { par[v][i + 1] = par[par[v][i]][i]; } } } int lca(int v, int u) { if (dep[v] < dep[u]) swap(v, u); int gap = dep[v] - dep[u]; for (int i = l - 1; i >= 0; i--) { if ((1 << i) & gap) v = par[v][i]; } if (v == u) return v; for (int i = l - 1; i >= 0; i--) { int pv = par[v][i]; int pu = par[u][i]; if (pv != pu) v = pv, u = pu; } return par[v][0]; } int length(int v, int u) { int a = lca(v, u); return dep[v] + dep[u] - dep[a] * 2; } T dist(int v, int u) { int a = lca(v, u); return dists[v] + dists[u] - dists[a] * 2; } };
lの正当性
lが境界値でも正しい値になるかどうかは、2進数で考えると理解しやすいです。
頂点数\(N\)のライングラフ（一直線上に頂点が並び、隣り合う頂点の間に辺が張られたグラフ）を考えます。両端のどちらかの頂点を根としたとき、もう一方の端の頂点の深さは\(N - 1\)です。これが頂点数\(N\)のときの深さの最大値です。
\(N = 8\)のときの最大の深さは7であり、これは2進数表記で100 + 10 + 1と表せるためl = 3で十分です。一方\(N = 9\)のときは、最大の深さは8であり、l = 3で表せる整数の範囲を超えるため、lは最低でも4以上必要です。
\(N = 1\)のとき、l = 0なので少し不安になります。意味的には、親を持つ頂点が存在しないので正しく、実装的には、根に限り親頂点の初期化の対象外としているため正しいです。そしてlca()の最後のreturn par[v][0];が実行されることもありません。
機能を追加する際は、l = 0のときでも正しく動作するよう注意する必要があります。
使用例
最小共通祖先 | グラフライブラリ | Aizu Online Judge
int main() { int N; cin >> N; auto g = LowestCommonAncestor<>(N); for (int i = 0; i < N; i++) { int k; cin >> k; for (int j = 0; j < k; j++) { int c; cin >> c; g.add_edge(i, c); } } g.build(); int Q; cin >> Q; for (int i = 0; i < Q; i++) { int v, u; cin >> v >> u; cout << g.lca(v, u) << '\n'; } return 0; }
オイラーツアーによる実装
オイラーツアーというのは、頂点を訪問順に並べたものです。訪問順というのは、任意の頂点からDFSを行ったときの訪問順のことです。頂点と限定してしまいましたが、訪れた「辺」でもオイラーツアーと呼びます。
以下の実装では、vs配列がオイラーツアーを表します。DFSを一度行うだけでオイラーツアーを構築できるのは注目すべき点です。
struct LowestCommonAncestor { using P = pair<int, int>; const int INF = 1e9; int n; int m; int segn; vector<vector<int>> es; vector<int> id; vector<int> vs; vector<int> depth; vector<P> seg; LowestCommonAncestor(int n) : n(n), es(n), id(n), m(n * 2 - 1), vs(m), depth(m) {} void add_edge(int v, int u) { es[v].push_back(u); es[u].push_back(v); } void dfs(int v, int p, int d, int &k) { id[v] = k; vs[k] = v; depth[k++] = d; for (auto u : es[v]) { if (u == p) continue; dfs(u, v, d + 1, k); vs[k] = v; depth[k++] = d; } } void seginit() { segn = 1; while (segn < m) segn <<= 1; seg.assign(segn * 2, P(INF, INF)); for (int i = 0; i < m; i++) { seg[i + segn] = P(depth[i], i); } for (int i = segn - 1; i > 0; i--) { seg[i] = min(seg[i * 2 + 0], seg[i * 2 + 1]); } } void build() { int k = 0; dfs(0, -1, 0, k); seginit(); } P query(int a, int b) { P L(INF, INF), R(INF, INF); for (a += segn, b += segn; a < b; a >>= 1, b >>= 1) { if (a & 1) L = min(L, seg[a++]); if (b & 1) R = min(seg[--b], R); } return min(L, R); } int lca(int v, int u) { if (id[v] > id[u]) swap(v, u); return vs[query(id[v], id[u] + 1).second]; } };
セグメント木（以下セグ木）とオイラーツアーの実装がひとつの構造体にまとめて記述されているので、少し理解しづらいかもしれません。ということで、セグ木を外部化した実装も載せておきます。
ちなみにセグ木は、範囲に関する情報を高速で取得したいときに使われるデータ構造です。通常の発想では、範囲の大きさが\(N\)のとき、範囲の情報の取得には最低でも\(O(N)\)かかりそうですが、セグ木を使うと\(O(\log N)\)まで減らすことができます。
ここでいう「範囲」は、横に並べられたものに対して使っています。グラフの範囲とか木の範囲といった掴みどころのない範囲に対してセグ木は使えませんが、オイラーツアーのような、木の頂点を横に並べたものに対してはセグ木が使えます（「そのままでは扱えないものを扱える形にする」というのは、競プロではよくある話です）。
セグメント木部分を分離した実装
template <typename T> struct SegmentTree { using F = function<T(T, T)>; vector<T> v; F f; T e; int n; SegmentTree(int size, F f, T e) : f(f), e(e) { n = 1; while (n < size) n <<= 1; v.resize(n * 2, e); } void set(int k, T x) { v[k + n] = x; } void build() { for (int i = n - 1; i > 0; i--) { v[i] = f(v[i * 2 + 0], v[i * 2 + 1]); } } void update(int k, T x) { v[k += n] = x; while (k >>= 1) v[k] = f(v[k * 2 + 0], v[k * 2 + 1]); } T query(int a, int b) { T l = e, r = e; for (a += n, b += n; a < b; a >>= 1, b >>= 1) { if (a & 1) l = f(l, v[a++]); if (b & 1) r = f(v[--b], r); } return f(l, r); } }; struct LowestCommonAncestor { using P = pair<int, int>; const int INF = 1e9; int n; int m; vector<vector<int>> es; vector<int> id; vector<int> vs; vector<int> depth; SegmentTree<P> seg; LowestCommonAncestor(int n) : n(n), es(n), id(n), m(n * 2 - 1), vs(m), depth(m), seg(m, [](P a, P b) { return min(a, b); }, P(INF, INF)) {} void add_edge(int v, int u) { es[v].push_back(u); es[u].push_back(v); } void dfs(int v, int p, int d, int &k) { id[v] = k; vs[k] = v; depth[k++] = d; for (auto u : es[v]) { if (u == p) continue; dfs(u, v, d + 1, k); vs[k] = v; depth[k++] = d; } } void build() { int k = 0; dfs(0, -1, 0, k); for (int i = 0; i < m; i++) seg.set(i, P(depth[i], i)); seg.build(); } int lca(int v, int u) { if (id[v] > id[u]) swap(v, u); return vs[seg.query(id[v], id[u] + 1).second]; } };
pair<int, int>に対してもmin()が使えるのは初めて知りました。
使用例
最小共通祖先 | グラフライブラリ | Aizu Online Judge
int main() { int N; cin >> N; auto g = LowestCommonAncestor(N); for (int i = 0; i < N; i++) { int k; cin >> k; for (int j = 0; j < k; j++) { int c; cin >> c; g.add_edge(i, c); } } g.build(); int Q; cin >> Q; for (int i = 0; i < Q; i++) { int v, u; cin >> v >> u; cout << g.lca(v, u) << '\n'; } return 0; }
HL分解による実装
HL分解は、木という扱いにくいものを扱いやすくするための手段のひとつです。HL分解により、木がパスの集合になります。パスは分岐がなく一次元のデータとして扱えるため、セグ木に載せることが可能です（今回この性質は利用しません）。
今回の目的はLCAを高速に求めることですが、このときに利用する性質は「HL分解後の木の高さは高々\(O(\log N)\)」です。
struct HeavyLightDecomposition { int n; vector<vector<int>> es; vector<int> par; vector<int> dep; vector<int> head; vector<int> next; vector<int> vid; HeavyLightDecomposition(int n) : n(n), es(n), par(n), dep(n), head(n), next(n, -1), vid(n) {} void add_edge(int v, int u) { es[v].push_back(u); es[u].push_back(v); } int dfs(int v = 0, int d = 0, int p = -1) { if (p != -1) par[v] = p; dep[v] = d; int mx = 0; int cnt = 1; for (int u : es[v]) { if (u == p) continue; int c = dfs(u, d + 1, v); cnt += c; if (mx < c) mx = c, next[v] = u; } return cnt; } void bfs(int v = 0) { int k = 0; queue<int> q; q.emplace(v); while (!q.empty()) { int h = q.front(); q.pop(); for (int v = h; v != -1; v = next[v]) { vid[v] = k++; head[v] = h; for (int u : es[v]) { if (u == par[v] || u == next[v]) continue; q.push(u); } } } } void build() { dfs(); bfs(); } int lca(int v, int u) { while (true) { if (vid[u] > vid[v]) swap(u, v); if (head[u] == head[v]) return u; v = par[head[v]]; } } };
「分解」といっても、別に辺を消去しているわけではないことがわかります。「グループ分け」をイメージしたほうがわかりやすいかもしれません。head[v]は、頂点\(v\)が属するグループの番号が入ります。UnionFindのpar[v]みたいなものです。
head[v]はもうひとつの意味を持ち、それは「グループ内で最も根に近い頂点」という意味です。したがって、par[head[v]]は「隣接している別グループ」を表します。HL分解後の木の高さは高々\(O(\log N)\)であるため、lca()で行われている「隣接しているグループに移動し続ける処理」の回数は\(O(\log N)\)に抑えられます。これがLCAを高速に求められる理由です。
vector<int>型の変数が5つありますが、これらはすべて頂点に関する情報を持つ変数です。以下の実装とほぼ同じと考えるとイメージしやすいかもしれません。
struct Vertex { int vid, dep; Vertex *par, *head, *next; }; vector<Vertex> vs;
使用例
最小共通祖先 | グラフライブラリ | Aizu Online Judge
int main() { int N; cin >> N; auto g = HeavyLightDecomposition(N); for (int i = 0; i < N; i++) { int k; cin >> k; for (int j = 0; j < k; j++) { int c; cin >> c; g.add_edge(i, c); } } g.build(); int Q; cin >> Q; for (int i = 0; i < Q; i++) { int v, u; cin >> v >> u; cout << g.lca(v, u) << '\n'; } return 0; }
Link-Cut Treeによる実装
struct LinkCutTree { struct Node { Node *l, *r, *p; int k; Node(int k) : k(k), l(nullptr), r(nullptr), p(nullptr) {} bool is_root() { return !p || (p->l != this && p->r != this); } }; int n; vector<Node *> nodes; LinkCutTree(int n) : n(n), nodes(n) { for (int i = 0; i < n; i++) { nodes[i] = new Node(i); } } void rotr(Node *t) { auto *x = t->p, *y = x->p; if ((x->l = t->r)) t->r->p = x; t->r = x, x->p = t; if ((t->p = y)) { if (y->l == x) y->l = t; if (y->r == x) y->r = t; } } void rotl(Node *t) { auto *x = t->p, *y = x->p; if ((x->r = t->l)) t->l->p = x; t->l = x, x->p = t; if ((t->p = y)) { if (y->l == x) y->l = t; if (y->r == x) y->r = t; } } void splay(Node *t) { while (!t->is_root()) { auto *p = t->p; if (p->is_root()) { if (p->l == t) rotr(t); else rotl(t); continue; } auto *q = p->p; if (q->l == p) { if (p->l == t) rotr(p), rotr(t); else rotl(t), rotr(t); } else { if (p->r == t) rotl(p), rotl(t); else rotr(t), rotl(t); } } } int expose(int k) { Node *t = nodes[k]; Node *rp = nullptr; for (auto *cur = t; cur; cur = cur->p) { splay(cur); cur->r = rp; rp = cur; } splay(t); return rp->k; } void link(int v, int u) { expose(v); expose(u); nodes[v]->p = nodes[u]; nodes[u]->r = nodes[v]; } int lca(int v, int u) { expose(v); return expose(u); } };
使用例
最小共通祖先 | グラフライブラリ | Aizu Online Judge
int main() { int N; cin >> N; LinkCutTree lct(N); for (int i = 0; i < N; i++) { int k; cin >> k; for (int j = 0; j < k; j++) { int u; cin >> u; lct.link(u, i); } } int Q; cin >> Q; for (int i = 0; i < Q; i++) { int v, u; cin >> v >> u; cout << lct.lca(v, u) << '\n'; } return 0; }
まとめ
LCAを求める方法を4つ紹介しました。
- ダブリング：各頂点が\(2 ^ i\)個上にある親の情報を持つことで、「2つの頂点\(v\)と\(u\)が同じになるまで上に辿り続ける」処理を高速化しています。
- オイラーツアー：木をオイラーツアーに変換し、LCAを求める問題を範囲の最小値を求める問題に言い換えています。これはセグ木を使って高速に求められます。
- HL分解：HL分解後の木の高さは\(O(\log N)\)であることを利用して高速化します。
- Link-Cut Tree：よくわかりませんが\(O(\log N)\)になります。
計算量はすべて\(O(\log N)\)です。方法が違うのに計算量は一致するのが面白いです。
