yukicoder No.274 The Wall 解説
問題
2-SATについての前提
2-SAT問題の解き方
- 論理和を論理包含に変換
- 1を元にimplication graphを作る
- 2のimplication graphを強連結成分分解する
- 3の結果、ある要素Xとその否定￢Xが強連結成分であるなら、その2-SATを満たす解は存在しない
式変換について
論理包含P→Qとは、Pが偽またはQが真であることを意味する。論理包含を利用することで論理和は、以下の様に書き換えることが可能。
P∨Q ⇔ ￢P→Q∧￢Q→P
implication graphについて
論理包含をの関係をグラフにおこしたものをimplication graphと言う。
implication graphにおいて強連結である頂点の真偽値は全て等しいという性質がある。（これは参考に載せた、2-SATを解いた時のメモに解説がある）
強連結である頂点は互いに行き来可能なので、(P, Q)が強連結であるなら、P→Q∧Q→P、つまりPとQは等しい。
この性質から、もしある変数Xとその否定￢Xが同じ強連結な集合に属するなら、そのimplication graphは矛盾することになり、その2-SATには解が無い事になる。
解説
1つの区間を2-SATにおける1つの変数と捉える。また、その変数の否定とは、区間を180度回転させたものとする。
すると、ある2つの区間が重なる可能性があるなら、その2つの区間を指す変数をA, Bとして、A∨Bという論理和の関係が成り立つ。（どちらか一方しか取れないため）
全ての2頂点間について、重なる可能性があるかを見ていく。（2頂点間の関係は180度回転させる場合を加味すると4パターンあることに注意する。）
上記の論理和の関係を、論理包含に変換しimplication graphを構築する。後は前提に記述した2-SATの解き方で解ける。
解答コード
struct StronglyConnectedComponents { int n; vector<vector<int>> G, rG; vector<int> order, component; vector<bool> used; void dfs(int v) { used[v] = 1; for (auto nv : G[v]) { if (!used[nv]) dfs(nv); } order.push_back(v); } void rdfs(int v, int k) { component[v] = k; for (auto nv : rG[v]) { if (component[nv] < 0) rdfs(nv, k); } } StronglyConnectedComponents(vector<vector<int>> &_G) { n = _G.size(); G = _G; rG.resize(n); component.assign(n, -1); used.resize(n); for (int v = 0; v < n; v++) { for (auto nv : G[v]) rG[nv].push_back(v); } for (int v = 0; v < n; v++) if (!used[v]) dfs(v); int k = 0; reverse(order.begin(), order.end()); for (auto v : order) if (component[v] == -1) rdfs(v, k), k++; } bool is_same(int u, int v) { return component[u] == component[v]; } vector<vector<int>> rebuild() { int N = *max_element(component.begin(), component.end()) + 1; vector<vector<int>> rebuildedG(N); set<pair<int, int>> connected; for (int v = 0; v < N; v++) { for (auto nv : G[v]) { if (component[v] != component[nv] and !connected.count(pair(v, nv))) { connected.insert(pair(v, nv)); rebuildedG[component[v]].push_back(component[nv]); } } } return rebuildedG; } }; #define rep(i, n) for (long long i = (long long)(0); i < (long long)(n); ++i) #define irep(i, m, n) for (long long i = (long long)(m); i < (long long)(n); ++i) using ll = long long; template<class t> using vc=vector<t>; signed main() { cin.tie( 0 ); ios::sync_with_stdio( false ); int n,m; cin>>n>>m; vc<int> l(n), r(n); rep(i,n) cin>>l[i]>>r[i]; vc<vc<int>> G(2*n); rep(i,n){ irep(j,i+1,n){ int l1=l[i], r1=r[i], l2=l[j], r2=r[j]; int L1=m-r1-1, R1=m-l1-1, L2=m-r2-1, R2=m-l2-1; if(l1<=l2 and l2<=r1 or l1<=r2 and r2<=r1 or l2<=l1 and l1<=r2){ // A∨B <=> ￢A=>B∧￢B=>A G[i+n].push_back(j); G[j+n].push_back(i); } if(l1<=L2 and L2<=r1 or l1<=R2 and R2<=r1 or L2<=l1 and l1<=R2){ // A∨￢B <=> ￢A=>￢B∧B=>A G[i+n].push_back(j+n); G[j].push_back(i); } if(L1<=l2 and l2<=R1 or L1<=r2 and r2<=R1 or l2<=L1 and L1<=r2){ // ￢A∨B <=> A=>B∧￢B=>￢A G[i].push_back(j); G[j+n].push_back(i+n); } if(L1<=L2 and L2<=R1 or L1<=R2 and R2<=R1 or L2<=L1 and L1<=R2){ // ￢A∨￢B <=> A=>￢B∧B=>￢A G[i].push_back(j+n); G[j].push_back(i+n); } } } auto scc = StronglyConnectedComponents(G); rep(i,n){ if(scc.is_same(i,i+n)){ cout<<"NO"<<endl; return 0; } } cout<<"YES"<<endl; }
参考
論理包含に関する解説
強連結成分分解の実装と個人的なメモ
前提
有向グラフGの頂点の部分集合をSとする。
- Sが強連結であるとは、Sの任意の2頂点間が行き来可能であること
- Sが強連結成分であるとは、Sに他のどの頂点を追加してもそれ以上強連結にならないこと
強連結成分分解とは、強連結成分を1つの頂点にまとめること。こうしてまとめられた頂点でグラフを再構築すると、その有向グラフはDAGになる。
- グラフをDAGに再構築し、DAGに対してDPするのは典型
- 強連結成分分解を使って解くことができる有名問題として2-SATがある
強連結成分分解のアルゴリズムは単純で2回DFSを行うだけ。
1回目のDFSで各頂点に帰りがけの順を記録する。
2回目のDFSは元のグラフの辺を逆辺に置き換えたものに対して行う。1回目で採番した帰りがけの遅い順にDFSする。既に訪れた頂点を避けながらDFSすることで、その頂点から行き来できる頂点集合が1つの強連結成分となる。
2回のDFSで何をしているのか
これは参考に載せた、目指せグラフマスターのスライドP.11によくまとめられている。
- 強連結な集合は辺を逆辺に置き換えても強連結のまま
- DFSの帰りがけ順が遅い => 完成したDAGのトポロジカル順序において早い
- 逆辺のDFSなら今見ている頂点から（2のトポロジカル順序における）後の強連結成分を訪れない
3が重要で、このトポロジカル順序を無視して逆辺のグラフにDFSすると、ある連結成分から、自身よりトポロジカル順序が前の連結成分へ移動してしまう可能性があり、本来連結成分ではない集合を連結成分としてしまう可能性がある。
実装
struct StronglyConnectedComponents { int n; vector<vector<int>> G, rG; vector<int> order, component; vector<bool> used; void dfs(int v) { used[v] = 1; for (auto nv : G[v]) { if (!used[nv]) dfs(nv); } order.push_back(v); } void rdfs(int v, int k) { component[v] = k; for (auto nv : rG[v]) { if (component[nv] < 0) rdfs(nv, k); } } StronglyConnectedComponents(vector<vector<int>> &_G) { n = _G.size(); G = _G; rG.resize(n); component.assign(n, -1); used.resize(n); for (int v = 0; v < n; v++) { for (auto nv : G[v]) rG[nv].push_back(v); } for (int v = 0; v < n; v++) if (!used[v]) dfs(v); int k = 0; reverse(order.begin(), order.end()); for (auto v : order) if (component[v] == -1) rdfs(v, k), k++; } /// 頂点(u, v)が同じ強連結成分に含まれるか bool is_same(int u, int v) { return component[u] == component[v]; } /// 強連結成分を1つのノードに潰したグラフを再構築する vector<vector<int>> rebuild() { int N = *max_element(component.begin(), component.end()) + 1; vector<vector<int>> rebuildedG(N); set<pair<int, int>> connected; for (int v = 0; v < N; v++) { for (auto nv : G[v]) { if (component[v] != component[nv] and !connected.count(pair(v, nv))) { connected.insert(pair(v, nv)); rebuildedG[component[v]].push_back(component[nv]); } } } return rebuildedG; } };
検証用問題
参考
p.82から強連結成分分解の解説
https://hcpc-hokudai.github.io/archive/graph_scc_001.pdf
www.slideshare.net
ダブリングによる接尾辞配列（Suffix-Array）の実装と個人的なメモ
AC Libraryに、より高速なライブラリが存在する。
以下に載せた実装は、自分の理解を深めるために実装したもの。（蟻本のものにアレンジを加えclass化した。計算量は変わらず、Manber & Myers の O(n log2 n) 時間のアルゴリズム）
前提
接尾辞配列（Suffix-Array）とは、ある文字列の接尾辞（Suffix）を辞書順にソートしたもの。
ナイーブにSuffix-Arrayを構築しようとすると、文字列をsとして、空間計算量 |s|^2、時間計算 |s|^2 log |s| となる。後述するが、Suffix-Arrayを構築するのは文字列中に含まれる部分文字列の検索を高速に行うためであり、ナイーブなSuffix-Arrayの構築の時間計算量は、ナイーブな検索の時間計算量(|s|^2)を上回ってしまうので、このままでは構築する意味があまりない。また空間計算量も、扱う文字列の長さが105以上の問題ではMLEになってしまう。
そのため、Suffix-Array構築時の計算量を改善する必要がある。
Suffix-Arrayでもつ値は部分文字列ではなくSuffixの始点のindexにする。(これで空間計算量は|s|になる)
各indexからのSuffixの辞書順の順位が求められれば、その順位を使ってSuffix-Arrayのindexをソートすることができるので、rank配列を作り、これをダブリングによって段階的に更新していく。
以下は、rank配列を更新する実装部分についての簡単な説明。
rank配列の初期化
rankはそのindexの文字コードによって初期化している。rank[n](空文字)だけは、-1で初期化。
for (int i = 0; i <= n; i++) { suffix_array[i] = i; rank[i] = i < n ? (int)s[i] : -1; }
確定した2*k文字まで長さのrankを元にSuffix-Arrayをソートする
k=1から始めて、各indexから2*k文字まで見た時の順位をrank配列に格納していく。
nを対象の文字列の長さとして、2*k>=nまで更新を行えば全てのSuffixのrankを求めたことになる。
ループ回数が log n 回、各シーケンス内のソートが n log n 時間なので、時間計算量は n log2 n
vector<int> tmp(n + 1); // 一度tmp配列に for (k = 1; k <= n; k *= 2) { sort(suffix_array.begin(), suffix_array.end(), compare); tmp[suffix_array[0]] = 0; for (int i = 1; i <= n; i++) { tmp[suffix_array[i]] = tmp[suffix_array[i - 1]] + (compare(suffix_array[i - 1], suffix_array[i]) ? 1 : 0); } for (int i = 0; i <= n; i++) { rank[i] = tmp[i]; } }
ソートに使う比較関数
k=1から始めて、2*kまでの長さの部分文字列の大小比較を行う。
auto compare = [&](int a, int b) { if (rank[a] != rank[b]) return rank[a] < rank[b]; // 前半k文字のrankが異なるなら、後半k文字を比較するまでもなく大小が確定する else { // 前半k文字のrankが同等なら、後半のk文字のrankの大小で2*k文字の大小が決まる // +k文字目が、文字列をはみ出すなら空文字としてrank=-1とする auto rank_ak = a + k <= n ? rank[a + k] : -1; auto rank_bk = b + k <= n ? rank[b + k] : -1; return rank_ak < rank_bk; } };
実装
class SuffixArray { private: string s; int n, k; vector<int> suffix_array; vector<int> rank; void create_suffix_array() { for (int i = 0; i <= n; i++) { suffix_array[i] = i; rank[i] = i < n ? (int)s[i] : -1; } auto compare = [&](int a, int b) { if (rank[a] != rank[b]) return rank[a] < rank[b]; else { auto rank_ak = a + k <= n ? rank[a + k] : -1; auto rank_bk = b + k <= n ? rank[b + k] : -1; return rank_ak < rank_bk; } }; vector<int> tmp(n + 1); for (k = 1; k <= n; k *= 2) { sort(suffix_array.begin(), suffix_array.end(), compare); tmp[suffix_array[0]] = 0; for (int i = 1; i <= n; i++) { tmp[suffix_array[i]] = tmp[suffix_array[i - 1]] + (compare(suffix_array[i - 1], suffix_array[i]) ? 1 : 0); } for (int i = 0; i <= n; i++) { rank[i] = tmp[i]; } } } public: SuffixArray(const string &_s) { s = _s; n = _s.size(); suffix_array.resize(n + 1); rank.resize(n + 1); create_suffix_array(); } // 文字列sの部分文字列として、パラメータで渡した文字列tが現れるかを返す bool contains(const string &t) { int l = 0, r = n; while (r - l > 1) { int mid = (l + r) / 2; int index = suffix_array[mid]; if (s.compare(index, t.size(), t) < 0) l = mid; else r = mid; } return s.compare(suffix_array[r], t.size(), t) == 0; } };
Suffix-Arrayを使った文字列検索
上記の実装のcontainsという関数は、文字列sの中にtが現れるかを判定している。
Suffix-Arrayがあることで二分探索によって部分文字列の検索を行えるので、O(|t| log |s|) 時間で判定を行える。
検証用問題
Suffix-Arrayを使った高速文字列検索によって解くことができる。（上記のcontainsにより検証可能）
ABC168 F - . (Single Dot) 解説
制約から座標圧縮をすべきであることは想像しやすいが、この問題はかなり難しいと思った。
前提
- 座標圧縮
解説
座標圧縮を行う問題でよく見かけるパターンは以下の2つである。
- いくつかの長方形が与えられて、長方形が重なる領域の面積の総和はいくつか？
- いくつかの線分が与えられて、線分で区切られた領域の数はいくつか？
今回の問題は、線分が与えられてそれにより区切られた領域の面積を求めるというもので、上記の2点とは微妙に異なる。
1に類する問題は、座標圧縮後の状態をグリッドにマッピングすると、各マスと面積が対応するので対象のマスの面積をそれぞれ求めてその総和を求めればよい。
2に類する問題は、マスと線分を対応させてグリッドにマッピングし、BFSすることで求めることができる。（この時、グリッドの1マスは面積に対応しない）
今回の問題では、マスを線分に対応させつつ、移動できるマスからその面積を求める必要がある。
線分をグリッドにマッピングしたとき、線分間に空間があるかをわかるようにするため、与えられる座標を全て2倍にする。
実際にサンプル1で試すとわかるが、座標圧縮しなくても、与えられる線分の情報をそのままグリッドにマッピングすると、線分間に空間があるのかわからない。例えば、x軸方向に延びる線分で、x=0, x=1, x=2 の3つの線分をグリッドにマッピングすると、当然3行を連続して塗りつぶしてしまい、この行間に空間が存在するのかがわからない。
そのため、各座標を2倍することを考える。x=0, x=1, x=2 → x=0, x=2, x=4 このようにすると、各線分間には必ず1以上の差が生まれるため、適切に座標圧縮する（線分の座標の値+1を空白として登録する。前提参照）ことでグリッドに線分間の空間をマッピングすることができる。
グリッドのマスから面積を求める方法については結論から書くと、マス (x, y) が空白でかつ x, y が共に奇数であるマスだけを数えるようにすると良い。
あるマスが空白であるなら、そのマスの面積は (隣のY座標 - 今のY座標) * (隣のX座標 - 今のX座標) で求めることができる。
これは、前工程で線分の各座標値を2倍にし、その隣の座標を空白としたため、(x, y) のいずれかが偶数のマスは線分間の空白をマッピングするためのバッファであり、そのマスの面積は線分間の面積に対応しないため。
実装
#define rep(i, n) for (long long i = (long long)(0); i < (long long)(n); ++i) #define all(v) v.begin(), v.end() using ll = long long; template<class t> using vc=vector<t>; template <typename T> vector<T> compress(vector<T> &v1, vector<T> &v2, vector<T> &v3) { vector<T> vals; int n = v1.size(), m = v3.size(); for (int i = 0; i < n; i++) { vals.push_back(v1[i]); vals.push_back(v2[i]); vals.push_back(v1[i] + 1); vals.push_back(v2[i] + 1); } for (int i = 0; i < m; i++) { vals.push_back(v3[i]); vals.push_back(v3[i] + 1); } vals.push_back(0); vals.push_back(1); sort(vals.begin(), vals.end()); vals.erase(unique(vals.begin(), vals.end()), vals.end()); for (int i = 0; i < n; i++) { v1[i] = lower_bound(vals.begin(), vals.end(), v1[i]) - vals.begin(); v2[i] = lower_bound(vals.begin(), vals.end(), v2[i]) - vals.begin(); } for (int i = 0; i < m; i++) { v3[i] = lower_bound(vals.begin(), vals.end(), v3[i]) - vals.begin(); } return vals; } void inf() { cout << "INF" << endl; exit(0); } const ll dy[] = {0, 1, 0, -1, -1, 1, 1, -1}; const ll dx[] = {1, 0, -1, 0, 1, 1, -1, -1}; inline bool inside(ll y, ll x, ll H, ll W) { return (y >= 0 && x >= 0 && y < H && x < W); } signed main() { cin.tie( 0 ); ios::sync_with_stdio( false ); int n, m; cin >> n >> m; vc<ll> x1(n), x2(n), y1(n), x3(m), y2(m), y3(m); rep(i, n) cin >> x1[i] >> x2[i] >> y1[i]; rep(i, m) cin >> x3[i] >> y2[i] >> y3[i]; // 全ての座標の値を2倍 rep(i, n) { x1[i] *= 2; x2[i] *= 2; y1[i] *= 2; } rep(i, m) { x3[i] *= 2; y2[i] *= 2; y3[i] *= 2; } // 座標圧縮 auto X = compress(x1, x2, x3); auto Y = compress(y2, y3, y1); int h = X.size(), w = Y.size(); // X軸方向、Y軸方向、別々にimos法を行う vc<vc<int>> g1(h, vc<int>(w)); // X軸方向 vc<vc<int>> g2(h, vc<int>(w)); // Y軸方向 rep(i, n) { g1[x1[i]][y1[i]]++; g1[x2[i] + 1][y1[i]]--; } rep(j, w) { rep(i, h - 1) { g1[i + 1][j] += g1[i][j]; } } rep(i, m) { g2[x3[i]][y2[i]]++; g2[x3[i]][y3[i] + 1]--; } rep(i, h) { rep(j, w - 1) { g2[i][j + 1] += g2[i][j]; } } // BFS で移動可能な領域の面積の総和を求める。h * w の範囲外に出られるなら INF vc<vc<bool>> visited(h, vc<bool>(w)); queue<pair<int, int>> q; int sx = lower_bound(X.begin(), X.end(), 0) - X.begin(); int sy = lower_bound(Y.begin(), Y.end(), 0) - Y.begin(); q.emplace(sx, sy); visited[sx][sy] = 1; ll ans = 0; while (!q.empty()) { int x, y; tie(x, y) = q.front(); q.pop(); if (x % 2 == 1 and y % 2 == 1) ans += (X[x + 1] / 2 - (X[x] - 1) / 2) * (Y[y + 1] / 2 - (Y[y] - 1) / 2); rep(i, 4) { int nx = x + dx[i], ny = y + dy[i]; if (inside(nx, ny, h, w) and g1[nx][ny] == 0 and g2[nx][ny] == 0 and visited[nx][ny] == 0) { visited[nx][ny] = 1; q.emplace(nx, ny); } else if (!inside(nx, ny, h, w)) { inf(); } } } cout << ans << endl; }
座標圧縮に関する解説と練習問題
前提
一次元座標圧縮はBITによる転倒数の数え上げでお馴染み。
二次元の場合は、x軸、y軸を独立に捉えることでほぼ同様に求めることができる。
二次元座標圧縮における注意点は、登場する座標の1つ隣の座標も登録する点にある。これは登場する座標だけで圧縮を行うと、元のグリッドの線分（や長方形）間に存在した空白を潰してしまう可能性があるため。（この図解が参考1のスライドにある）
座標圧縮の基本的な実装
- パラメータに渡した座標をmutableに扱い、圧縮後の座標に置き換える
- 座標圧縮が必要になる問題は、線分の始点と終点や長方形（3Dなら直方体）の対角など、同数の2つの座標が与えられる形式が多く、そのため引数に2つのリストを取るようにすると便利なことが多い（練習問題での compress の呼び出し元を参照）
- 戻り値として、引数に渡した座標（とその1つ隣の座標）の圧縮前の座標をソートし重複を排除したリスト
valsを返す
上記の様にすることで、圧縮後の値xから圧縮前の値を復元したい場合、vals[x] と記述するだけで簡単に行うことができる。
実際の使い方は以下の練習問題の解法を参照。
template<typename T> vector<T> compress(vector<T> &v1, vector<T> &v2){ vector<T> vals; int n = v1.size(); for(int i=0; i<n; i++){ for(int j=0; j<=1; j++){ vals.push_back(v1[i]+j); vals.push_back(v2[i]+j); } } sort(vals.begin(), vals.end()); vals.erase(unique(vals.begin(), vals.end()), vals.end()); for(int i=0; i<n; i++){ v1[i] = lower_bound(vals.begin(), vals.end(), v1[i]) - vals.begin(); v2[i] = lower_bound(vals.begin(), vals.end(), v2[i]) - vals.begin(); } return vals; }
練習問題
座標圧縮が必要になる問題では、座標圧縮を行い、imos法により圧縮後の座標をグリッドにマッピングするという処理を書くことが多い。（そのため、ここでは多次元imos法の書き方を知っていることを前提とする）
そして以下の様なことを問う問題が多い。
- 座標平面と長方形が与えられて、長方形が重なる領域の面積の総和はいくつか？（三次元の問題もある）
- 座標平面と線分が与えられて、線分で区切られた領域の数はいくつか？
二次元座標圧縮
https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=DSL_4_A&lang=ja
問題形式1の問題で、圧縮した座標をグリッドにマッピングし、そのグリッドの1マス1マスの面積を求め、その総和を求めることで解を求めている。
signed main() { int n; cin>>n; vector<long long> x1(n), y1(n), x2(n), y2(n); for(int i=0; i<n; i++) cin>>x1[i]>>y1[i]>>x2[i]>>y2[i]; vector<long long> X = compress(x1, x2); vector<long long> Y = compress(y1, y2); int h = Y.size(), w = X.size(); vector<vector<int>> G(h, vector<int>(w)); // 二次元imos for(int i=0; i<n; i++){ G[y1[i]][x1[i]]++; G[y2[i]][x2[i]]++; G[y1[i]][x2[i]]--; G[y2[i]][x1[i]]--; } for(int i=0; i<h-1; i++){ for(int j=0; j<w; j++){ G[i+1][j] += G[i][j]; } } for(int i=0; i<h; i++){ for(int j=0; j<w-1; j++){ G[i][j+1] += G[i][j]; } } // グリッドにおける各セルの面積の総和を求めることで解を得る long long ans=0; for(int i=0; i<h; i++){ for(int j=0; j<w; j++){ if(G[i][j]>0){ ans += (Y[i+1]-Y[i]) * (X[j+1]-X[j]); } } } cout<<ans<<endl; }
問題形式2の問題。座標圧縮したグリッドに対してBFSを行い、その回数が解となる。
この問題では圧縮する座標にグリッドの両端を加えるため、引数に head, tail を追加している。
#define rep(i, n) for (long long i = (long long)(0); i < (long long)(n); ++i) template<class t> using vc=vector<t>; template <typename T> vector<T> compress(vector<T> &v1, vector<T> &v2, T head, T tail) { vector<T> vals; int n = v1.size(); for (int i = 0; i < n; i++) { vals.push_back(v1[i]); vals.push_back(v2[i]); if(v2[i] + 1 <= tail) vals.push_back(v2[i] + 1); } vals.push_back(head); vals.push_back(tail); sort(vals.begin(), vals.end()); vals.erase(unique(vals.begin(), vals.end()), vals.end()); for (int i = 0; i < n; i++) { v1[i] = lower_bound(vals.begin(), vals.end(), v1[i]) - vals.begin(); v2[i] = lower_bound(vals.begin(), vals.end(), v2[i]) - vals.begin(); } return vals; } const int dy[] = {0, 1, 0, -1}; const int dx[] = {1, 0, -1, 0}; inline bool inside(int y, int x, int H, int W) { return (y >= 0 && x >= 0 && y < H && x < W); } signed main() { cin.tie( 0 ); ios::sync_with_stdio( false ); int H,W; cin>>H>>W; int n; cin>>n; vc<int> x1(n),y1(n),x2(n),y2(n); rep(i,n) cin>>x1[i]>>y1[i]>>x2[i]>>y2[i]; auto X=compress(x1,x2,0,H); auto Y=compress(y1,y2,0,W); int h=X.size(), w=Y.size(); vc<vc<int>> G(h, vc<int>(w)); // 二次元imos rep(i,n){ G[x1[i]][y1[i]]++; G[x1[i]][y2[i]]--; G[x2[i]][y1[i]]--; G[x2[i]][y2[i]]++; } rep(i,h){ rep(j,w-1){ G[i][j+1]+=G[i][j]; } } rep(j,w){ rep(i,h-1){ G[i+1][j]+=G[i][j]; } } h--, w--; // BFS int ans=0; rep(i,h){ rep(j,w){ if(G[i][j]==0){ ans++; queue<pair<int,int>> q; q.emplace(i,j); G[i][j]=1; while(!q.empty()){ int x,y; tie(x,y)=q.front(); q.pop(); rep(k,4){ int nx=x+dx[k], ny=y+dy[k]; if(inside(nx,ny,h,w) and G[nx][ny]==0){ G[nx][ny]=1; q.emplace(nx,ny); } } } } } } cout<<ans<<endl; }
三次元座標圧縮
これも問題形式1の問題で三次元に拡張したもの。（三次元imos法の書き方の練習にもなる）
#define rep(i, n) for (long long i = (long long)(0); i < (long long)(n); ++i) using ll = long long; template<class t> using vc=vector<t>; signed main() { cin.tie( 0 ); ios::sync_with_stdio( false ); int n,K; cin>>n>>K; vc<ll> x1(n), y1(n), z1(n), x2(n), y2(n), z2(n); rep(i,n){ cin>>x1[i]>>y1[i]>>z1[i]>>x2[i]>>y2[i]>>z2[i]; } vc<ll> X=compress(x1,x2); vc<ll> Y=compress(y1,y2); vc<ll> Z=compress(z1,z2); int h=X.size(), w=Y.size(), d=Z.size(); vc<vc<vc<ll>>> G(h, vc<vc<ll>>(w, vc<ll>(d))); // 三次元imos rep(i,n){ G[x1[i]][y1[i]][z1[i]]++; G[x1[i]][y2[i]][z1[i]]--; G[x1[i]][y1[i]][z2[i]]--; G[x2[i]][y1[i]][z1[i]]--; G[x1[i]][y2[i]][z2[i]]++; G[x2[i]][y1[i]][z2[i]]++; G[x2[i]][y2[i]][z1[i]]++; G[x2[i]][y2[i]][z2[i]]--; } rep(i,h){ rep(j,w){ rep(k,d-1){ G[i][j][k+1]+=G[i][j][k]; } } } rep(i,h){ rep(k,d){ rep(j,w-1){ G[i][j+1][k]+=G[i][j][k]; } } } rep(k,d){ rep(j,w){ rep(i,h-1){ G[i+1][j][k]+=G[i][j][k]; } } } ll ans=0; rep(i,h){ rep(j,w){ rep(k,d){ if(G[i][j][k]>=K){ ans+=(X[i+1]-X[i])*(Y[j+1]-Y[j])*(Z[k+1]-Z[k]); } } } } cout<<ans<<endl; }
参考
参考1）
www.slideshare.net
参考2）実装を参考にさせていただきました。
Dijkstra法による最小費用流を求める実装（Primal-Dual）とポテンシャルに関する個人的なメモ
前提
基本的な考えは、Bellman-Ford法による最小費用流の求め方と同じである。
以下の実装では高速化のためBellman-Ford法の実装個所をDijkstra法に置き換えており、同時に負の辺が存在するグラフに対してDijkstra法を適用するため、ポテンシャルを導入している。
計算量: O(flow|E|log|V|)（最大流flow、頂点数V、辺数Eとする）
ポテンシャルについて
※ 以下の説明は、niuezさんのブログ（参考1）を参考に自分の理解を整理するために書いたものです。
最短距離の数式化
単一始点最短経路を求め済みのグラフにおいて、頂点iの最短距離をd[i]、頂点(i,j)の辺のコストをdist(i,j)とすると、辺が存在する頂点(i, j)間には以下のことが成り立つ。（これが最短距離を数式として落とし込んだものになる）
- d[s] = 0 とする。（sは始点)
- d[i] + dist(i, j) >= d[j] （iが始点に近く、jが終点に近い頂点とする）
ポテンシャルを導入し、グラフ（辺のコスト）を再構築する
各頂点にポテンシャルhを決めて、辺の距離を dist'(i, j) = dist(i, j) + h[i] - h[j] に置き換えたグラフG'を考える。 G'における各頂点の最短距離をd'とする。G'にも、最短距離の数式が成り立ち、更に2行目の様に式変換ができる。
d'[i] + dist'(i, j) >= d'[j]
⇔ d'[i] + h[i] + dist(i, j) >= d'[j] + h[j]
2行目から、d[i] = d'[i] + h[i] と見れば、d[i]は最短距離の式の形をしており、d[s] = d'[s] + h[s] より各頂点iについて d'[i] + h[i] - h[s] を計算すれば元の最短距離が求められるとわかる。
各辺を非負とするポテンシャルの求め方
ここまでで、ポテンシャルを導入して元のグラフの最短距離を求められることが分かった。次に、Dijkstra法を適用できるようにするためdist'を非負にしたいので、ポテンシャルに適切な値を設定することを考える。
dist(i, j) + h[i] - h[j] >=0
⇔ h[i] + dist(i, j) >= h[j]
2行目から最短距離の式の形をしているため、始点からの最短距離がdist'を非負にするためのポテンシャルであるとわかる。つまり1度Bellman-Ford法等により、各頂点の最短距離を求めておけば、以降その値をポテンシャルとすることでDijkstra法によっても最短距離を求めることができるようになる。
まとめ
- 一度Bellman-Ford法等によって最短距離を求められればポテンシャルも求まる
- ポテンシャルから非負の辺のみからなるグラフG'を構築でき、これに対してDijkstra法を適用できる
- 2で求めたG'における各頂点の最短距離d'から始点sのポテンシャルを減算すれば元のグラフにおける最短距離を求められる
ちなみに以下の最小費用流を求める実装においては、Dijkstra法によって1度目の最短距離を求めている。これは、初期状態のグラフにおいて移動可能な経路上に負の辺が存在しないことを前提としているため。
参考
参考1）ポテンシャルについて
参考2）負閉路が無い事とポテンシャルが存在することは同値であることの説明
検証用問題
実装
template <typename capacity_t, typename cost_t> class PrimalDual { private: const cost_t INF = 1e9; struct Edge { int to, rev; capacity_t cap; cost_t cost; Edge(int to, int rev, capacity_t cap, cost_t cost) : to(to), rev(rev), cap(cap), cost(cost) {} }; vector<vector<Edge>> G; vector<int> prev_v, prev_e; Edge &get_rev(Edge &edge) { return G[edge.to][edge.rev]; } public: PrimalDual(int n) { G.resize(n); prev_v.resize(n); prev_e.resize(n); } void add_edge(int from, int to, capacity_t cap, cost_t cost) { G[from].push_back(Edge(to, (int)G[to].size(), cap, cost)); G[to].push_back(Edge(from, (int)G[from].size() - 1, 0, -cost)); } cost_t get_min_cost_flow(int s, int t, capacity_t flow) { int n = G.size(); vector<cost_t> h(n); // ポテンシャル cost_t res = 0; while (flow > 0) { vector<cost_t> dist(n, INF); dist[s] = 0; priority_queue<pair<cost_t, int>, vector<pair<cost_t, int>>, greater<pair<cost_t, int>>> que; que.emplace(0, s); while (!que.empty()) { cost_t c; int v; tie(c, v) = que.top(); que.pop(); if (c > dist[v]) continue; for (int i = 0; i < G[v].size(); i++) { auto &edge = G[v][i]; if (edge.cap > 0 and dist[edge.to] > dist[v] + edge.cost + h[v] - h[edge.to]) { dist[edge.to] = dist[v] + edge.cost + h[v] - h[edge.to]; prev_v[edge.to] = v; prev_e[edge.to] = i; que.emplace(dist[edge.to], edge.to); } } } if (dist[t] == INF) return -1; for (int v = 0; v < n; v++) h[v] += dist[v]; capacity_t d = flow; for (int v = t; v != s; v = prev_v[v]) { d = min(d, G[prev_v[v]][prev_e[v]].cap); } flow -= d; res += d * h[t]; for (int v = t; v != s; v = prev_v[v]) { G[prev_v[v]][prev_e[v]].cap -= d; get_rev(G[prev_v[v]][prev_e[v]]).cap += d; } } return res; } };
Bellman-Ford法による最小費用流を求める実装
前提
最小費用流問題のグラフには、最大流問題で与えられる辺の情報に加えて辺のコストが付与される。
以下の実装では、このコストに対してS-T最短路を1つ見つけて目いっぱい流す、その経路が限界になったら次の最短路を見つけて目いっぱい流して、という事を求める流量に達するまで反復して行うことで最小費用流を求めている。
最小費用流問題を求める際にも残余グラフを構築するが、この際、逆辺のコストは元のコストを-1倍した値を設定することになる。
故に残余グラフには負の辺を含むことになるため、そのままではDijkstra法を使うことができないため、最短路を求めるためにBellman-Ford法を利用している。（ポテンシャルを導入することで最短路を求める工程をDijkstra法に置き換える方法もある => 続き）
計算量: O(flow|V||E|)（最大流flow、頂点数V、辺数Eとする）
参考
P.48 以降参照
http://hos.ac/slides/20150319_flow.pdf
検証用問題
実装
template <typename capacity_t, typename cost_t> class MinCostFlow { private: const cost_t INF = 1e9; struct Edge { int to, rev; capacity_t cap; cost_t cost; Edge(int to, int rev, capacity_t cap, cost_t cost) : to(to), rev(rev), cap(cap), cost(cost) {} }; vector<vector<Edge>> G; vector<int> prev_v, prev_e; Edge& get_rev(Edge &edge) { return G[edge.to][edge.rev]; } public: MinCostFlow(int n) { G.resize(n); prev_v.resize(n); prev_e.resize(n); } void add_edge(int from, int to, capacity_t cap, cost_t cost) { G[from].push_back(Edge(to, (int)G[to].size(), cap, cost)); G[to].push_back(Edge(from, (int)G[from].size() - 1, 0, -cost)); } cost_t get_min_cost_flow(int s, int t, capacity_t flow) { int n = G.size(); cost_t res = 0; while (flow > 0) { vector<cost_t> dist(n, INF); dist[s] = 0; bool update = true; while (update) { update = false; for (int v = 0; v < n; v++) { if (dist[v] == INF) continue; for (int i = 0; i < G[v].size(); i++) { auto &edge = G[v][i]; if (edge.cap > 0 and dist[edge.to] > dist[v] + edge.cost) { dist[edge.to] = dist[v] + edge.cost; prev_v[edge.to] = v; prev_e[edge.to] = i; update = true; } } } } if (dist[t] == INF) return -1; capacity_t d = flow; for (int v = t; v != s; v = prev_v[v]) { d = min(d, G[prev_v[v]][prev_e[v]].cap); } flow -= d; res += d * dist[t]; for (int v = t; v != s; v = prev_v[v]) { G[prev_v[v]][prev_e[v]].cap -= d; get_rev(G[prev_v[v]][prev_e[v]]).cap += d; } } return res; } };
