// 森林火災の伐採問題: 参加者指定のビームサーチ // usage: main [--th 4] [--width 10] [--seed 0] < input > output // --th : 乾燥度がこの値以下のマスは常に湿気、それより大きいマスは常に乾燥とみなす // --width : dp[伐採本数] ごとに保持する状態数(ビーム幅) #include using namespace std; constexpr int N = 50, NN = N * N; int K; int TH = 3, WIDTH = 50; uint64_t SEED = 0; int d[NN]; bool dryBase[NN]; // 伐採前に乾燥とみなすマス mt19937_64 rng; struct State { bitset cut; int cnt = 0; long long eval = 0; // Σ (連結成分サイズ)^2 }; struct UF { int p[NN], sz[NN]; void init() { for (int i = 0; i < NN; i++) p[i] = i, sz[i] = 1; } 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; inline bool isDry(const State& s, int c) { return dryBase[c] && !s.cut[c]; } // 乾燥マスの 4 近傍連結成分を union-find で求める void buildUF(const State& s) { uf.init(); for (int c = 0; c < NN; c++) { if (!isDry(s, c)) continue; int i = c / N, j = c % N; if (j + 1 < N && isDry(s, c + 1)) uf.unite(c, c + 1); if (i + 1 < N && isDry(s, c + N)) uf.unite(c, c + N); } } long long evaluate(const State& s) { buildUF(s); long long e = 0; for (int c = 0; c < NN; c++) if (isDry(s, c) && uf.find(c) == c) e += (long long)uf.sz[c] * uf.sz[c]; return e; } // 乾燥マス上の 4 近傍 BFS 距離(到達不能は -1)。最大距離を返す int dist[NN], que[NN]; int bfs(const State& s, int src) { fill(dist, dist + NN, -1); int h = 0, t = 0, mx = 0; dist[src] = 0, que[t++] = src; while (h < t) { int c = que[h++], i = c / N, j = c % N; mx = dist[c]; int nb[4] = {i > 0 ? c - N : -1, i + 1 < N ? c + N : -1, j > 0 ? c - 1 : -1, j + 1 < N ? c + 1 : -1}; for (int u : nb) if (u >= 0 && dist[u] < 0 && isDry(s, u)) dist[u] = dist[c] + 1, que[t++] = u; } return mx; } int pickFarthest(const State& s, int src) { int mx = bfs(s, src); vector cand; for (int c = 0; c < NN; c++) if (dist[c] == mx) cand.push_back(c); return cand[rng() % cand.size()]; } vector> dp; void insert(State&& st) { auto& b = dp[st.cnt]; for (auto& o : b) if (o.cut == st.cut) return; // 伐採集合が重複したら追加しない if ((int)b.size() < WIDTH) { b.push_back(std::move(st)); return; } auto worst = max_element(b.begin(), b.end(), [](const State& a, const State& c) { return a.eval < c.eval; }); if (st.eval < worst->eval) *worst = std::move(st); } // 最大成分 L 内で ban されたマスを除いた連結成分の Σ サイズ^2(BFS) int mark[NN], stampId = 0; long long ssWithin(const State& s, 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++], i = c / N, j = c % N; int nb[4] = {i > 0 ? c - N : -1, i + 1 < N ? c + N : -1, j > 0 ? c - 1 : -1, j + 1 < N ? c + 1 : -1}; for (int u : nb) if (u >= 0 && mark[u] != stampId && isDry(s, u)) mark[u] = stampId, que[t++] = u; } ss += (long long)t * t; } return ss; } // 端点からの距離がちょうど x の層を伐採したときの L 内の Σ サイズ^2 を全 x について求める。 // 距離 < x の部分は連結で、距離 > x の部分とは隣接しないので、 // 距離 > x の部分の成分を x の降順に union-find で育てて求める(evaluate と同じ値になる)。 UF uf2; bool inR[NN]; vector ssLayers(const vector& L, int diam) { int maxd = 0; for (int c : L) maxd = max(maxd, dist[c]); vector> byDist(maxd + 1); for (int c : L) byDist[dist[c]].push_back(c), inR[c] = false, uf2.p[c] = c, uf2.sz[c] = 1; long long ssGt = 0; auto addCell = [&](int c) { inR[c] = true; ssGt += 1; int i = c / N, j = c % N; int nb[4] = {i > 0 ? c - N : -1, i + 1 < N ? c + N : -1, j > 0 ? c - 1 : -1, j + 1 < N ? c + 1 : -1}; for (int u : nb) { if (u < 0 || !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 x = maxd; x > diam; x--) for (int c : byDist[x]) addCell(c); vector res(diam + 1); long long less = 0; vector cntLess(diam + 1); for (int x = 0; x <= diam; x++) cntLess[x] = less, less += byDist[x].size(); for (int x = diam; x >= 0; x--) { res[x] = ssGt + cntLess[x] * cntLess[x]; for (int c : byDist[x]) addCell(c); } for (int c : L) inR[c] = false; return res; } void expand(const State& s) { buildUF(s); int best = 0; vector roots; for (int c = 0; c < NN; c++) { if (!isDry(s, 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()]; // 最大成分(同率は一様ランダム)の代表元 vector L; for (int c = 0; c < NN; c++) if (isDry(s, c) && uf.find(c) == v) L.push_back(c); long long baseEval = s.eval - (long long)L.size() * L.size(); int prev = v; for (int r = 0; r < 5; r++) prev = v, v = pickFarthest(s, v); int ends[2] = {prev, v}; // 直径 s, t bfs(s, ends[0]); int diam = dist[ends[1]]; for (int e : ends) { bfs(s, e); vector ss = ssLayers(L, diam); vector> layer(diam + 1); for (int c : L) if (dist[c] <= diam) layer[dist[c]].push_back(c); for (int x = 0; x <= diam; x++) { vector add = layer[x]; if (add.empty()) continue; State ch = s; if (s.cnt + (int)add.size() > K) { shuffle(add.begin(), add.end(), rng); add.resize(K - s.cnt); for (int c : add) ch.cut[c] = 1; ch.eval = baseEval + ssWithin(ch, L, add); } else { for (int c : add) ch.cut[c] = 1; ch.eval = baseEval + ss[x]; } ch.cnt += add.size(); #ifdef CHECK_EVAL assert(ch.eval == evaluate(ch)); #endif insert(std::move(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); int n; cin >> n >> K; for (int i = 0; i < N; i++) { string row; cin >> row; for (int j = 0; j < N; j++) d[i * N + j] = row[j] - '0'; } for (int c = 0; c < NN; c++) dryBase[c] = d[c] > TH; dp.assign(K + 1, {}); State init; init.eval = evaluate(init); dp[0].push_back(init); for (int c = 0; c < K; c++) { vector cur = dp[c]; for (auto& s : cur) expand(s); } if (dp[K].empty()) { fprintf(stderr, "dp[K] is empty\n"); return 1; } auto& res = *min_element(dp[K].begin(), dp[K].end(), [](const State& a, const State& b) { return a.eval < b.eval; }); fprintf(stderr, "eval = %lld\n", res.eval); printf("%d\n", res.cnt); for (int c = 0; c < NN; c++) if (res.cut[c]) printf("%d %d\n", c / N, c % N); }