#include #include using namespace std; struct ForestRng { uint64_t x; explicit ForestRng(uint64_t s = 0x6a09e667f3bcc909ULL) : x(s) {} uint64_t operator()() { x ^= x << 7; x ^= x >> 9; return x; } int integer(int n) { return int((*this)() % uint64_t(n)); } double real() { return double((*this)() >> 11) * 0x1.0p-53; } }; struct ForestSamples { int N, V, maxSamples; vector d, vertices; vector> adj; vector degree, dry; vector order, low, sub, sepSum, sepSq, component; const uint8_t *dryNow = nullptr; const uint8_t *cutNow = nullptr; int timer = 0; ForestSamples( const vector& grid, int sampleCount = 512, uint64_t seed = 6523487656ULL ) : N((int)grid.size()), V(N * N), maxSamples(sampleCount), d(V), adj(V), degree(V), dry(size_t(V) * sampleCount), order(V), low(V), sub(V), sepSum(V), sepSq(V), component(V) { for (int i = 0; i < N; i++) { for (int j = 0; j < N; j++) { int u = i * N + j; d[u] = grid[i][j] - '0'; if (d[u]) vertices.push_back(u); } } for (int u : vertices) { int i = u / N, j = u % N; if (i && d[u - N]) adj[u][degree[u]++] = u - N; if (j && d[u - 1]) adj[u][degree[u]++] = u - 1; if (j + 1 < N && d[u + 1]) adj[u][degree[u]++] = u + 1; if (i + 1 < N && d[u + N]) adj[u][degree[u]++] = u + N; } // 各20サンプルで,各頂点を正確に2*d回乾燥させる. ForestRng rng(seed); const int block = 20; array idx; for (int base = 0; base < sampleCount; base += block) { int len = min(block, sampleCount - base); for (int u : vertices) { iota(idx.begin(), idx.end(), 0); int quota = d[u] * len / 10; if (d[u] * len % 10 && rng.integer(10) < d[u] * len % 10) { ++quota; } for (int t = 0; t < quota; t++) { swap(idx[t], idx[t + rng.integer(len - t)]); dry[size_t(base + idx[t]) * V + u] = 1; } } } } void dfs(int u, int parent) { order[u] = low[u] = ++timer; component[timer - 1] = u; int sz = 1, ss = 0, sq = 0; for (int k = 0; k < degree[u]; k++) { int v = adj[u][k]; if (v == parent || cutNow[v] || !dryNow[v]) continue; if (!order[v]) { dfs(v, u); sz += sub[v]; low[u] = min(low[u], low[v]); if (low[v] >= order[u]) { ss += sub[v]; sq += sub[v] * sub[v]; } } else { low[u] = min(low[u], order[v]); } } sub[u] = sz; sepSum[u] = ss; sepSq[u] = sq; } // gainsを指定すると,全頂点の削除利得も計算する. int64_t evaluate( const vector& cuts, int samples, vector* gains = nullptr, int offset = 0 ) { assert(samples + offset <= maxSamples); cutNow = cuts.data(); if (gains) fill(gains->begin(), gains->end(), 0); int64_t total = 0; for (int s = offset; s < offset + samples; s++) { dryNow = dry.data() + size_t(s) * V; fill(order.begin(), order.end(), 0); timer = 0; if (gains) { for (int u : vertices) { if (!dryNow[u] || cutNow[u] || order[u]) continue; int begin = timer; dfs(u, -1); int sz = timer - begin; total += sz * sz; for (int z = begin; z < timer; z++) { int v = component[z]; int rem = sz - 1 - sepSum[v]; (*gains)[v] += int64_t(sz) * sz - sepSq[v] - rem * rem; } } } else { // 利得不要なら,通常の連結成分探索だけを行う. for (int u : vertices) { if (!dryNow[u] || cutNow[u] || order[u]) continue; int head = timer, begin = timer; component[timer++] = u; order[u] = 1; while (head < timer) { int v = component[head++]; for (int k = 0; k < degree[v]; k++) { int w = adj[v][k]; if (!order[w] && dryNow[w] && !cutNow[w]) { order[w] = 1; component[timer++] = w; } } } int sz = timer - begin; total += sz * sz; } } } return total; } }; // 全て伐採した状態から,木を復元する構築法. struct ReverseBuilder { ForestSamples& m; ForestRng rng; int S, V, offset; vector parent, size; vector cut; ReverseBuilder( ForestSamples& model, int samples, uint64_t seed, int off = 0 ) : m(model), rng(seed), S(samples), V(model.V), offset(off), parent(samples * V, -1), size(samples * V), cut(V) { for (int u : m.vertices) cut[u] = 1; } int root(int s, int u) { int off = s * V; while (parent[off + u] != u) { parent[off + u] = parent[off + parent[off + u]]; u = parent[off + u]; } return u; } // u自身は常に乾燥すると仮定して周囲との合併コストを求め, // 最後に乾燥度を掛ける.共通係数1/(10*S)は省略する. double cost(int u) { int64_t total = 0; for (int s = 0; s < S; s++) { const uint8_t* dry = m.dry.data() + size_t(s + offset) * V; int roots[4], count = 0; int merged = 1, oldSquares = 0; for (int k = 0; k < m.degree[u]; k++) { int v = m.adj[u][k]; if (cut[v] || !dry[v]) continue; int r = root(s, v); bool duplicate = false; for (int j = 0; j < count; j++) { if (roots[j] == r) duplicate = true; } if (duplicate) continue; roots[count++] = r; int z = size[s * V + r]; merged += z; oldSquares += z * z; } total += merged * merged - oldSquares; } return double(total) * m.d[u]; } void restore(int u) { cut[u] = 0; for (int s = 0; s < S; s++) { const uint8_t* dry = m.dry.data() + size_t(s + offset) * V; if (!dry[u]) continue; int off = s * V; parent[off + u] = u; size[off + u] = 1; for (int k = 0; k < m.degree[u]; k++) { int v = m.adj[u][k]; if (cut[v] || !dry[v]) continue; int a = root(s, u), b = root(s, v); if (a == b) continue; if (size[off + a] < size[off + b]) swap(a, b); parent[off + b] = a; size[off + a] += size[off + b]; } } } vector run(int K) { vector remaining = m.vertices; while ((int)remaining.size() > K) { int chosen = -1; double best = 1e100; bool exhaustive = remaining.size() <= 300; int count = exhaustive ? (int)remaining.size() : 40; for (int k = 0; k < count; k++) { int j = exhaustive ? k : rng.integer((int)remaining.size()); double value = cost(remaining[j]) * (0.999 + 0.002 * rng.real()); if (value < best) { best = value; chosen = j; } } restore(remaining[chosen]); remaining[chosen] = remaining.back(); remaining.pop_back(); } return cut; } }; struct Solver { using Clock = chrono::steady_clock; Clock::time_point start; double limit; ForestSamples model; ForestRng rng; int K; vector> archive; Solver(const vector& grid, int k, double seconds) : start(Clock::now()), limit(seconds), model(grid, 1800), rng(987654321), K(min(k, (int)model.vertices.size())) {} double elapsed() const { return chrono::duration(Clock::now() - start).count(); } int64_t validation(const vector& cut, int samples = 200) { return model.evaluate(cut, samples, nullptr, 1000); } vector refine(vector cut, double deadline) { vector removed; for (int u : model.vertices) { if (cut[u]) removed.push_back(u); } vector gain(model.V); int S = 20; int64_t score = model.evaluate(cut, S); int64_t bestScore = score; auto best = cut; double begin = elapsed(); double duration = max(1e-6, deadline - begin); while (elapsed() < deadline && !removed.empty()) { double progress = (elapsed() - begin) / duration; int target = progress < .22 ? 20 : progress < .52 ? 40 : progress < .78 ? 80 : 160; if (target != S) { archive.push_back(best); cut = best; removed.clear(); for (int u : model.vertices) { if (cut[u]) removed.push_back(u); } S = target; score = bestScore = model.evaluate(cut, S); } int index = rng.integer((int)removed.size()); int add = removed[index]; cut[add] = 0; int64_t restored = model.evaluate(cut, S, &gain); int del = -1; for (int u : model.vertices) { if (cut[u] || u == add) continue; if (del < 0 || gain[u] > gain[del] || (gain[u] == gain[del] && (rng() & 7) == 0)) { del = u; } } if (del < 0) { cut[add] = 1; break; } int64_t next = restored - gain[del]; double phase = fmod(progress * 4, 1.0); double temperature = score * .0015 * pow(.02, phase); if (progress > .88) temperature *= .05; if (next <= score || rng.real() < exp(double(score - next) / max(1e-9, temperature))) { cut[del] = 1; removed[index] = del; score = next; if (score < bestScore) { bestScore = score; best = cut; } } else { cut[add] = 1; } } archive.push_back(best); return best; } vector solve() { if (K == (int)model.vertices.size()) { vector cut(model.V); for (int u : model.vertices) cut[u] = 1; return cut; } if (K == 0) return vector(model.V); vector>> candidates; int attempts = 0; do { ReverseBuilder builder( model, 20, 987654321ULL + 55434ULL * attempts, 0 ); auto cut = builder.run(K); candidates.emplace_back(validation(cut), move(cut)); ++attempts; } while (elapsed() < limit * .25); sort(candidates.begin(), candidates.end()); auto best = candidates[0].second; archive.push_back(best); best = refine(best, limit * .91); vector> ranked; for (int i = 0; i < (int)archive.size(); i++) { ranked.emplace_back(validation(archive[i]), i); } sort(ranked.begin(), ranked.end()); int64_t bestScore = LLONG_MAX; for (int j = 0; j < min(3, (int)ranked.size()); j++) { int idx = ranked[j].second; int64_t value = validation(archive[idx], 800); if (value < bestScore) { bestScore = value; best = archive[idx]; } } cerr << "starts=" << attempts << " elapsed=" << elapsed() << " seconds\n"; return best; } }; int main() { ios::sync_with_stdio(false); cin.tie(nullptr); int N, K; if (!(cin >> N >> K)) return 0; vector grid(N); for (auto& row : grid) cin >> row; double seconds = 1.85; if (const char* s = getenv("TIME_LIMIT")) { seconds = max(.1, atof(s)); } Solver solver(grid, K, seconds); auto cut = solver.solve(); cout << accumulate(cut.begin(), cut.end(), 0) << '\n'; for (int u = 0; u < N * N; u++) { if (cut[u]) { cout << u / N << ' ' << u % N << '\n'; } } }