// 森林火災の伐採問題: 参加者指定のビームサーチ(差分更新版) // usage: main [--th 4] [--width 10] [--seed 0] < input > output // --th : 乾燥度がこの値以下のマスは常に湿気、それより大きいマスは常に乾燥とみなす // --width : dp[伐採本数] ごとに保持する状態数(ビーム幅) // // 状態はコピーせず、探索木(各ノードは親と「追加で伐採したマス」の差分だけを持つ)で管理する。 // dp[c] を展開するときは、dp[c] の状態の祖先だけからなる部分木を DFS し、 // 1 つの盤面に差分を適用・巻き戻ししながら各状態を展開する。 // 探索木は伐採本数の飛び(世代スキップ)がある遷移でもそのまま扱える。 // // 盤面は周囲 1 マスの番兵(常に乾燥していない)を付けた W x W で持つ。 // マス (i, j) の番号は (i + 1) * W + (j + 1)。行優先の走査順は元の盤面と同じ。 #pragma GCC optimize("O3") #include using namespace std; constexpr int N = 50, W = N + 2, WW = W * W; constexpr int DIRS[4] = {-W, W, -1, 1}; int K; int TH = 3, WIDTH = 100; uint64_t SEED = 0; int d[WW]; bool dryBase[WW]; // 伐採前に乾燥とみなすマス(番兵は false) uint64_t zob[WW]; // 伐採集合の重複判定用ハッシュ mt19937_64 rng; vector cells; // 盤内のマス(行優先) // 現在の盤面(差分を適用・巻き戻しして動かす) bool dryNow[WW]; // dryBase かつ伐採されていない int curCnt = 0; long long curEval = 0; uint64_t curHash = 0; inline bool isDry(int c) { return dryNow[c]; } inline void setCut(int c, bool v) { dryNow[c] = dryBase[c] && !v; } struct UF { int p[WW], sz[WW]; int find(int x) { while (p[x] != x) x = p[x] = p[p[x]]; return x; } void unite(int a, int b) { a = find(a), b = find(b); if (a == b) return; if (sz[a] < sz[b]) swap(a, b); p[b] = a, sz[a] += sz[b]; } } uf, uf2; // 乾燥マスの 4 近傍連結成分を union-find で求める void buildUF() { for (int c : cells) uf.p[c] = c, uf.sz[c] = 1; for (int c : cells) { if (!isDry(c)) continue; if (isDry(c + 1)) uf.unite(c, c + 1); if (isDry(c + W)) uf.unite(c, c + W); } } long long evaluateFull() { buildUF(); long long e = 0; for (int c : cells) if (isDry(c) && uf.find(c) == c) e += (long long)uf.sz[c] * uf.sz[c]; return e; } // 乾燥マス上の 4 近傍 BFS 距離。distSt[c] != bfsId のマスは未到達。 // 到達したマスは que[0, qlen) に BFS 順に並ぶ。最大距離を返す int dist[WW], distSt[WW], que[WW], qlen = 0, bfsId = 0; int bfs(int src) { ++bfsId; int h = 0, t = 0; dist[src] = 0, distSt[src] = bfsId, que[t++] = src; while (h < t) { int c = que[h++], nd = dist[c] + 1; for (int dd : DIRS) { int u = c + dd; if (distSt[u] != bfsId && isDry(u)) dist[u] = nd, distSt[u] = bfsId, que[t++] = u; } } qlen = t; return dist[que[t - 1]]; } int pickFarthest(int src) { int mx = bfs(src); static vector cand; cand.clear(); for (int k = qlen - 1; k >= 0 && dist[que[k]] == mx; k--) cand.push_back(que[k]); sort(cand.begin(), cand.end()); // マス番号順(候補の並びを固定するため) return cand[rng() % cand.size()]; } // 最大成分 L 内で ban されたマスを除いた連結成分の Σ サイズ^2(BFS) int mark[WW], stampId = 0; long long ssWithin(const vector& L, const vector& ban) { ++stampId; for (int c : ban) mark[c] = stampId; // ban 済みと訪問済みを同じ印で表す long long ss = 0; for (int src : L) { if (mark[src] == stampId) continue; int h = 0, t = 0; mark[src] = stampId, que[t++] = src; while (h < t) { int c = que[h++]; for (int dd : DIRS) { int u = c + dd; if (mark[u] != stampId && isDry(u)) mark[u] = stampId, que[t++] = u; } } ss += (long long)t * t; } return ss; } // 距離ごとのマス一覧(L を距離でバケットソートしたもの)。 // 層 x のマスは byDist[layerBeg[x], layerBeg[x + 1]) に L の順で並ぶ int byDist[WW], layerBeg[WW + 2]; int bucketByDist(const vector& L) { int maxd = 0; for (int c : L) maxd = max(maxd, dist[c]); fill(layerBeg, layerBeg + maxd + 2, 0); for (int c : L) layerBeg[dist[c] + 1]++; for (int x = 0; x <= maxd; x++) layerBeg[x + 1] += layerBeg[x]; static int pos[WW + 1]; copy(layerBeg, layerBeg + maxd + 1, pos); for (int c : L) byDist[pos[dist[c]]++] = c; return maxd; } // 端点からの距離がちょうど x の層を伐採したときの L 内の Σ サイズ^2 を全 x について求める。 // 距離 < x の部分は連結で、距離 > x の部分とは隣接しないので、 // 距離 > x の部分の成分を x の降順に union-find で育てて求める(evaluateFull と同じ値になる)。 bool inR[WW]; long long ssRes[WW]; void ssLayers(const vector& L, int maxd, int diam) { for (int c : L) inR[c] = false, uf2.p[c] = c, uf2.sz[c] = 1; long long ssGt = 0; auto addCell = [&](int c) { inR[c] = true; ssGt += 1; for (int dd : DIRS) { int u = c + dd; if (!inR[u]) continue; int a = uf2.find(c), b = uf2.find(u); if (a == b) continue; ssGt += 2LL * uf2.sz[a] * uf2.sz[b]; uf2.unite(a, b); } }; for (int k = layerBeg[diam + 1]; k < layerBeg[maxd + 1]; k++) addCell(byDist[k]); for (int x = diam; x >= 0; x--) { long long less = layerBeg[x]; ssRes[x] = ssGt + less * less; for (int k = layerBeg[x]; k < layerBeg[x + 1]; k++) addCell(byDist[k]); } for (int c : L) inR[c] = false; } // 探索木のノード(展開された状態) struct Node { int parent; int cnt; long long eval; uint64_t hash; int diffBegin, diffLen; // diffPool[diffBegin, diffBegin + diffLen) が親から追加で伐採したマス }; vector nodes; vector diffPool; // まだ展開していない候補 struct Cand { int parent; long long eval; uint64_t hash; vector diff; }; // dp[c]: eval が大きい順のヒープ(先頭が最悪)+ 伐採集合ハッシュの集合 struct Bucket { vector heap; unordered_set hashes; }; vector dp; auto worseFirst = [](const Cand& a, const Cand& b) { return a.eval < b.eval; }; // dp[c] に入るかどうか(入らない候補は作らずに捨てる) inline bool admissible(int c, long long eval, uint64_t hash) { auto& b = dp[c]; if ((int)b.heap.size() >= WIDTH && eval >= b.heap.front().eval) return false; return !b.hashes.count(hash); // 伐採集合が重複したら追加しない } void insert(int c, Cand&& cd) { auto& b = dp[c]; if ((int)b.heap.size() < WIDTH) { b.hashes.insert(cd.hash); b.heap.push_back(std::move(cd)); push_heap(b.heap.begin(), b.heap.end(), worseFirst); return; } pop_heap(b.heap.begin(), b.heap.end(), worseFirst); b.hashes.erase(b.heap.back().hash); b.hashes.insert(cd.hash); b.heap.back() = std::move(cd); push_heap(b.heap.begin(), b.heap.end(), worseFirst); } void applyNode(int v) { const Node& n = nodes[v]; for (int k = 0; k < n.diffLen; k++) setCut(diffPool[n.diffBegin + k], true); curCnt = n.cnt, curEval = n.eval, curHash = n.hash; } void revertNode(int v) { const Node& n = nodes[v]; for (int k = 0; k < n.diffLen; k++) setCut(diffPool[n.diffBegin + k], false); const Node& p = nodes[n.parent]; curCnt = p.cnt, curEval = p.eval, curHash = p.hash; } // 現在の盤面(ノード self)を展開する void expand(int self) { buildUF(); int best = 0; static vector roots; roots.clear(); for (int c : cells) { if (!isDry(c) || uf.find(c) != c) continue; if (uf.sz[c] > best) best = uf.sz[c], roots.clear(); if (uf.sz[c] == best) roots.push_back(c); } if (roots.empty()) return; // 乾燥マスがない int v = roots[rng() % roots.size()]; // 最大成分(同率は一様ランダム)の代表元 static vector L; L.clear(); for (int c : cells) if (isDry(c) && uf.find(c) == v) L.push_back(c); long long baseEval = curEval - (long long)L.size() * L.size(); int prev = v; for (int r = 0; r < 5; r++) prev = v, v = pickFarthest(v); int ends[2] = {prev, v}; // 直径 s, t // 最後の pickFarthest は s から BFS しているので、dist は s からの距離になっている int diam = dist[ends[1]]; static vector add; for (int e : ends) { if (e == ends[1]) bfs(e); int maxd = bucketByDist(L); ssLayers(L, maxd, diam); for (int x = 0; x <= diam; x++) { int beg = layerBeg[x], len = layerBeg[x + 1] - beg; if (len == 0) continue; add.assign(byDist + beg, byDist + beg + len); long long eval; if (curCnt + len > K) { shuffle(add.begin(), add.end(), rng); add.resize(K - curCnt); eval = baseEval + ssWithin(L, add); } else { eval = baseEval + ssRes[x]; } uint64_t hash = curHash; for (int c : add) hash ^= zob[c]; int nc = curCnt + add.size(); #ifdef CHECK_EVAL for (int c : add) setCut(c, true); assert(eval == evaluateFull()); for (int c : add) setCut(c, false); #endif if (!admissible(nc, eval, hash)) continue; insert(nc, Cand{self, eval, hash, add}); } } } // dp[c] の候補をノード化し、その祖先からなる部分木を DFS して展開する vector> childList; vector visitStamp; vector isTarget; int curStep = 0; void dfs(int v) { if (isTarget[v]) { expand(v); return; } // 対象どうしは祖先関係にならない(伐採本数が同じ) for (int ch : childList[v]) { applyNode(ch); dfs(ch); revertNode(ch); } } int main(int argc, char** argv) { for (int a = 1; a + 1 < argc; a += 2) { string k = argv[a]; if (k == "--th") TH = atoi(argv[a + 1]); else if (k == "--width") WIDTH = atoi(argv[a + 1]); else if (k == "--seed") SEED = strtoull(argv[a + 1], nullptr, 10); } rng.seed(SEED); { mt19937_64 zr(0x9e3779b97f4a7c15ULL); for (int c = 0; c < WW; c++) zob[c] = zr(); } int n; cin >> n >> K; for (int i = 0; i < N; i++) { string row; cin >> row; for (int j = 0; j < N; j++) { int c = (i + 1) * W + (j + 1); d[c] = row[j] - '0'; cells.push_back(c); } } for (int c : cells) dryBase[c] = d[c] > TH, dryNow[c] = dryBase[c]; dp.assign(K + 1, {}); nodes.push_back({0, 0, evaluateFull(), 0, 0, 0}); // 根 curEval = nodes[0].eval; // 根を展開 expand(0); int finalBest = -1; for (int c = 1; c <= K; c++) { auto& b = dp[c]; if (b.heap.empty()) continue; // 候補をノード化 vector targets; for (auto& cd : b.heap) { Node nd{cd.parent, c, cd.eval, cd.hash, (int)diffPool.size(), (int)cd.diff.size()}; diffPool.insert(diffPool.end(), cd.diff.begin(), cd.diff.end()); targets.push_back(nodes.size()); nodes.push_back(nd); } vector().swap(b.heap); unordered_set().swap(b.hashes); if (c == K) { for (int t : targets) if (finalBest < 0 || nodes[t].eval < nodes[finalBest].eval) finalBest = t; break; } // 対象の祖先だけからなる部分木を作る int M = nodes.size(); childList.resize(M); visitStamp.resize(M, -1); isTarget.assign(M, 0); ++curStep; vector used; for (int t : targets) { isTarget[t] = 1; int v = t; while (v != 0 && visitStamp[v] != curStep) { visitStamp[v] = curStep; used.push_back(nodes[v].parent); childList[nodes[v].parent].push_back(v); v = nodes[v].parent; } } dfs(0); for (int u : used) childList[u].clear(); } if (finalBest < 0) { fprintf(stderr, "dp[K] is empty\n"); return 1; } // 根から finalBest までの差分を集める vector out; for (int v = finalBest; v != 0; v = nodes[v].parent) for (int k = 0; k < nodes[v].diffLen; k++) out.push_back(diffPool[nodes[v].diffBegin + k]); sort(out.begin(), out.end()); fprintf(stderr, "eval = %lld\n", nodes[finalBest].eval); printf("%d\n", (int)out.size()); for (int c : out) printf("%d %d\n", c / W - 1, c % W - 1); }