#pragma GCC optimize("O3") #include using namespace std; // tools.zip/src/lib.rs の採点処理を、先頭10サンプルについて移植する。 // 採用する操作は、指定された「伐採位置を1つ、別の乾燥木へ移す」のみ。 // スコアが厳密に改善した操作だけを採用する山登り法。 constexpr int SAMPLES = 10; constexpr double SEARCH_TIME_LIMIT_MS = 1900.0; constexpr int STALL_TRIALS = 10000; constexpr int FINAL_SAMPLES = 1000; constexpr uint64_t GAMMA = 0x9E3779B97F4A7C15ULL; uint64_t splitmix64(uint64_t &state) { state += GAMMA; uint64_t z = state; z = (z ^ (z >> 30)) * 0xBF58476D1CE4E5B9ULL; z = (z ^ (z >> 27)) * 0x94D049BB133111EBULL; return z ^ (z >> 31); } struct Board { int n, k, v; vector d; vector cut; vector> neighbors; vector degree; Board(int n_, int k_, vector d_) : n(n_), k(k_), v(n * n), d(move(d_)), cut(v, 0), neighbors(v), degree(v, 0) { for (int r = 0; r < n; ++r) { for (int c = 0; c < n; ++c) { const int x = r * n + c; if (r > 0) neighbors[x][degree[x]++] = x - n; if (r + 1 < n) neighbors[x][degree[x]++] = x + n; if (c > 0) neighbors[x][degree[x]++] = x - 1; if (c + 1 < n) neighbors[x][degree[x]++] = x + 1; } } } // Rust の sample_seed(input) と同一。 uint64_t sample_seed() const { uint64_t h = 1469598103934665603ULL; for (uint8_t x : d) { h ^= uint64_t('0' + x); h *= 1099511628211ULL; } return h ^ (uint64_t(n) * 1000003ULL + uint64_t(k) * 998244353ULL); } // Rust の sample_dry(input, sample_index) と同一。 vector generate_dry(uint64_t sample_index) const { uint64_t state = sample_seed() + uint64_t(sample_index) * uint64_t((v + 1) / 2) * GAMMA; vector dry(v, 0); for (int x = 0; x < v; x += 2) { const uint64_t r = splitmix64(state); const uint64_t lo_threshold = uint64_t(double(d[x]) * (4294967296.0 / 10.0)); dry[x] = uint8_t((r & 0xFFFFFFFFULL) < lo_threshold); if (x + 1 < v) { const uint64_t hi_threshold = uint64_t(double(d[x + 1]) * (4294967296.0 / 10.0)); dry[x + 1] = uint8_t((r >> 32) < hi_threshold); } } return dry; } }; // 伐採なしの1サンプルにおける、連結成分のサイズの二乗和。 long long score_without_cuts(const Board &b, const vector &dry) { vector seen(b.v, 0); vector queue; queue.reserve(b.v); long long total = 0; for (int x = 0; x < b.v; ++x) { if (!dry[x] || seen[x]) continue; queue.clear(); queue.push_back(x); seen[x] = 1; for (size_t h = 0; h < queue.size(); ++h) { int u = queue[h]; for (int j = 0; j < b.degree[u]; ++j) { int w = b.neighbors[u][j]; if (dry[w] && !seen[w]) { seen[w] = 1; queue.push_back(w); } } } total += 1LL * queue.size() * queue.size(); } return total; } // 現在の連結成分ごとに、各マスを取り除いたときの分割サイズを保持する。 // 候補の評価では連結成分の探索を行わず、採用された場合だけ再構築する。 // 評価関数と、交換候補の選択・採用条件は変更しない。 struct SampleComponents { const Board &b; vector dry; vector component, component_size; vector tin, tout, low, parent, subtree_size; vector rest_size; vector> separating_child; vector separating_count; vector deletion_delta; long long q = 0; int timer = 0; SampleComponents(const Board &board, vector dry_) : b(board), dry(move(dry_)), component(b.v, -1), tin(b.v, -1), tout(b.v), low(b.v), parent(b.v, -1), subtree_size(b.v), rest_size(b.v), separating_child(b.v), separating_count(b.v, 0), deletion_delta(b.v) { rebuild(); } void rebuild() { fill(tin.begin(), tin.end(), -1); fill(component.begin(), component.end(), -1); fill(parent.begin(), parent.end(), -1); component_size.clear(); timer = 0; q = 0; // DFSの木と到達可能な辺を使って、各頂点を取り除いた場合の分割を調べる。 auto dfs = [&](auto &&self, int u, int id) -> void { component[u] = id; tin[u] = low[u] = timer++; subtree_size[u] = 1; for (int j = 0; j < b.degree[u]; ++j) { int w = b.neighbors[u][j]; if (!dry[w] || b.cut[w] || w == parent[u]) continue; if (tin[w] < 0) { parent[w] = u; self(self, w, id); subtree_size[u] += subtree_size[w]; low[u] = min(low[u], low[w]); } else { low[u] = min(low[u], tin[w]); } } tout[u] = timer; }; for (int u = 0; u < b.v; ++u) { if (!dry[u] || b.cut[u] || tin[u] >= 0) continue; int id = (int)component_size.size(); dfs(dfs, u, id); int sz = subtree_size[u]; component_size.push_back(sz); q += 1LL * sz * sz; } for (int u = 0; u < b.v; ++u) { if (component[u] < 0) continue; const int sz = component_size[component[u]]; int apart = 0; long long after_sum = 0; int count = 0; for (int j = 0; j < b.degree[u]; ++j) { int w = b.neighbors[u][j]; if (parent[w] == u && low[w] >= tin[u]) { separating_child[u][count++] = w; apart += subtree_size[w]; after_sum += 1LL * subtree_size[w] * subtree_size[w]; } } separating_count[u] = (uint8_t)count; rest_size[u] = sz - 1 - apart; after_sum += 1LL * rest_size[u] * rest_size[u]; deletion_delta[u] = after_sum - 1LL * sz * sz; } } // メモリ上の成分情報だけで、交換による正確な評価値の差分を返す。 long long evaluate_swap(int add, int remove) const { const int old_id = component[remove]; // sample上で非乾燥なら -1 long long delta = old_id < 0 ? 0 : deletion_delta[remove]; if (!dry[add]) return delta; int used_keys[4]; int used_count = 0; long long joined_size = 1, before_join = 0; for (int j = 0; j < b.degree[add]; ++j) { int w = b.neighbors[add][j]; if (w == remove || !dry[w] || b.cut[w]) continue; int id = component[w]; int key = id, sz = component_size[id]; if (id == old_id) { // removeによって切り離されるDFSの子の部分木に含まれるかを調べる。 key = old_id; sz = rest_size[remove]; for (int t = 0; t < separating_count[remove]; ++t) { int child = separating_child[remove][t]; if (tin[child] <= tin[w] && tin[w] < tout[child]) { key = -1 - t; sz = subtree_size[child]; break; } } } bool unique = true; for (int t = 0; t < used_count; ++t) { if (used_keys[t] == key) { unique = false; break; } } if (!unique) continue; used_keys[used_count++] = key; joined_size += sz; before_join += 1LL * sz * sz; } return delta + joined_size * joined_size - before_join; } }; long long official_score_10(long long q0, long long q, int v) { // tools.zip の Simulation::e, e0, score と同じ計算。 const double e0 = (double)q0 / SAMPLES / v; const double e = (double)q / SAMPLES / v; return llround(1e6 * (1.0 + log2(e0 / max(e, 1e-12)))); } // 候補の比較に使う1000サンプルは、tools.zip と同じ生成方法で得る。 // 各行の乾燥マスを64bit整数に記録して、全候補について使い回す。 // 連結成分のサイズ二乗和を、各行内の連続する木をまとめて計算する。 struct FinalEvaluator { const Board &b; vector sample_rows; vector parent, size; explicit FinalEvaluator(const Board &board) : b(board), sample_rows((size_t)FINAL_SAMPLES * b.n), parent(b.v), size(b.v) { const uint64_t seed = b.sample_seed(); const uint64_t advance = uint64_t((b.v + 1) / 2) * GAMMA; vector thresholds(b.v); for (int x = 0; x < b.v; ++x) { thresholds[x] = uint64_t(double(b.d[x]) * (4294967296.0 / 10.0)); } if (b.n % 2 == 0) { // 偶数幅では64bit乱数の上位・下位32bitが同じ行に収まる。 for (int s = 0; s < FINAL_SAMPLES; ++s) { uint64_t state = seed + uint64_t(s) * advance; uint64_t *rows = sample_rows.data() + (size_t)s * b.n; for (int i = 0; i < b.n; ++i) { uint64_t mask = 0; const uint64_t *thr = thresholds.data() + i * b.n; for (int j = 0; j < b.n; j += 2) { const uint64_t rnd = splitmix64(state); mask |= uint64_t((rnd & 0xffffffffULL) < thr[j]) << j; mask |= uint64_t((rnd >> 32) < thr[j + 1]) << (j + 1); } rows[i] = mask; } } } else { // 奇数幅では乱数のペアが行境界をまたぐため、元と同じ順序で処理。 for (int s = 0; s < FINAL_SAMPLES; ++s) { uint64_t state = seed + uint64_t(s) * advance; uint64_t *rows = sample_rows.data() + (size_t)s * b.n; for (int x = 0; x < b.v; x += 2) { const uint64_t rnd = splitmix64(state); if ((rnd & 0xffffffffULL) < thresholds[x]) { rows[x / b.n] |= 1ULL << (x % b.n); } if (x + 1 < b.v && (rnd >> 32) < thresholds[x + 1]) { rows[(x + 1) / b.n] |= 1ULL << ((x + 1) % b.n); } } } } } inline int find(int x) { while (parent[x] != x) { parent[x] = parent[parent[x]]; x = parent[x]; } return x; } inline long long unite(int a, int c) { a = find(a); c = find(c); if (a == c) return 0; if (size[a] < size[c]) swap(a, c); long long delta = 2LL * size[a] * size[c]; parent[c] = a; size[a] += size[c]; return delta; } long long evaluate(const vector &cuts) { uint64_t cut_rows[64] = {}; for (int v : cuts) cut_rows[v / b.n] |= 1ULL << (v % b.n); long long total_q = 0; for (int s = 0; s < FINAL_SAMPLES; ++s) { const uint64_t *row = sample_rows.data() + (size_t)s * b.n; uint64_t prev_mask = 0, prev_starts = 0; long long q = 0; for (int i = 0; i < b.n; ++i) { const uint64_t mask = row[i] & ~cut_rows[i]; const uint64_t starts = mask & ~(mask << 1); uint64_t x = starts; while (x) { const int j = __builtin_ctzll(x); x &= x - 1; const uint64_t suffix = mask >> j; const int len = (~suffix) ? __builtin_ctzll(~suffix) : 64 - j; const int id = i * b.n + j; parent[id] = id; size[id] = len; q += 1LL * len * len; } uint64_t both_starts = (prev_mask & mask) & ~((prev_mask & mask) << 1); while (both_starts) { const int j = __builtin_ctzll(both_starts); both_starts &= both_starts - 1; const uint64_t prefix = ~0ULL >> (63 - j); const int prev_j = 63 - __builtin_clzll(prev_starts & prefix); const int now_j = 63 - __builtin_clzll(starts & prefix); q += unite((i - 1) * b.n + prev_j, i * b.n + now_j); } prev_mask = mask; prev_starts = starts; } total_q += q; } return total_q; } }; int main() { ios::sync_with_stdio(false); cin.tie(nullptr); int n, k; if (!(cin >> n >> k)) return 0; vector d(n * n); for (int i = 0; i < n; ++i) { string row; cin >> row; for (int j = 0; j < n; ++j) d[i * n + j] = uint8_t(row[j] - '0'); } Board board(n, k, move(d)); vector candidate_cells; candidate_cells.reserve(board.v); for (int x = 0; x < board.v; ++x) { if (board.d[x] >= 6) candidate_cells.push_back(x); } mt19937_64 rng((uint64_t)chrono::steady_clock::now().time_since_epoch().count()); shuffle(candidate_cells.begin(), candidate_cells.end(), rng); int m = min(k, (int)candidate_cells.size()); vector cut_cells(candidate_cells.begin(), candidate_cells.begin() + m); vector uncut_cells(candidate_cells.begin() + m, candidate_cells.end()); for (int x : cut_cells) board.cut[x] = 1; // 現在の伐採位置を引き継ぎ、各10サンプルの山登り終了時の最良解を保存。 vector> saved_cuts; const auto start = chrono::steady_clock::now(); #ifdef HC_STATS uint64_t trials = 0, accepted = 0, stalls = 0; #endif if (!cut_cells.empty() && !uncut_cells.empty()) { vector sample_states; sample_states.reserve(SAMPLES); long long q0_total = 0, q_total = 0; uint64_t next_sample_index = 0; auto load_samples = [&]() { sample_states.clear(); q0_total = q_total = 0; for (int s = 0; s < SAMPLES; ++s) { vector dry = board.generate_dry(next_sample_index++); q0_total += score_without_cuts(board, dry); sample_states.emplace_back(board, move(dry)); q_total += sample_states.back().q; } }; load_samples(); long long current_score = official_score_10(q0_total, q_total, board.v); int unsuccessful = 0; uniform_int_distribution cut_dist(0, (int)cut_cells.size() - 1); uniform_int_distribution uncut_dist(0, (int)uncut_cells.size() - 1); while (chrono::duration( chrono::steady_clock::now() - start).count() < SEARCH_TIME_LIMIT_MS) { const int ci = cut_dist(rng); const int ui = uncut_dist(rng); const int add = cut_cells[ci]; const int remove = uncut_cells[ui]; long long delta_total = 0; for (int s = 0; s < SAMPLES; ++s) { delta_total += sample_states[s].evaluate_swap(add, remove); } const long long candidate_score = official_score_10( q0_total, q_total + delta_total, board.v ); #ifdef HC_STATS ++trials; #endif if (candidate_score > current_score) { board.cut[add] = 0; board.cut[remove] = 1; for (int s = 0; s < SAMPLES; ++s) { if (sample_states[s].dry[add] || sample_states[s].dry[remove]) { sample_states[s].rebuild(); } } swap(cut_cells[ci], uncut_cells[ui]); q_total += delta_total; current_score = candidate_score; unsuccessful = 0; #ifdef HC_STATS ++accepted; #endif } else if (++unsuccessful == STALL_TRIALS) { saved_cuts.push_back(cut_cells); load_samples(); // 次の10サンプルを取り直し、現在の解から続ける。 current_score = official_score_10(q0_total, q_total, board.v); unsuccessful = 0; #ifdef HC_STATS ++stalls; #endif } } } saved_cuts.push_back(cut_cells); // 最後の10サンプルでの最良解も保存する。 // 1000サンプルをすべての保存解に共通で用い、期待燃焼数が最小の解を選ぶ。 // 同じ解が複数回保存された場合、評価結果を再利用する。 // 交換で配列の順序だけが変わる場合もあるので、各集合をソートして同一判定する。 for (auto &cuts : saved_cuts) sort(cuts.begin(), cuts.end()); sort(saved_cuts.begin(), saved_cuts.end()); saved_cuts.erase(unique(saved_cuts.begin(), saved_cuts.end()), saved_cuts.end()); FinalEvaluator evaluator(board); long long best_q = LLONG_MAX; vector best_cuts; for (const auto &cuts : saved_cuts) { long long q = evaluator.evaluate(cuts); if (q < best_q) { best_q = q; best_cuts = cuts; } } #ifdef HC_STATS cerr << "trials=" << trials << " accepted=" << accepted << " resamples=" << stalls << " saved=" << saved_cuts.size() << " final_q=" << best_q << " elapsed_ms=" << chrono::duration(chrono::steady_clock::now() - start).count() << '\n'; #endif cout << best_cuts.size() << '\n'; for (int x : best_cuts) cout << x / n << ' ' << x % n << '\n'; return 0; }