// https://github.com/nmnmnmnm04/ahc-template #include #include using namespace std; using ll = long long; using ld = long double; using u64 = unsigned long long; using u128 = unsigned __int128; #define rep(i, n) for (ll i = 0; i < (n); ++i) #define drep(i, n) for (ll i = (n) - 1; i >= 0; --i) #define nfor(i, s, n) for (ll i = (s); i < (n); ++i) #define fi first #define se second #define mp make_pair #define pb push_back #define endl "\n" #define nall(v) (v).begin(), (v).end() #define rall(v) (v).rbegin(), (v).rend() #define pcnt(x) __builtin_popcountll(x) #define next_p(v) next_permutation((v).begin(), (v).end()) template bool chmin(T& current, const U& candidate) { if (current <= candidate) return false; current = candidate; return true; } template bool chmax(T& current, const U& candidate) { if (current >= candidate) return false; current = candidate; return true; } const ll INF = 1LL << 61; // 先頭4方向は U, D, L, R、後ろ4方向は斜め const int di[8] = {-1, 1, 0, 0, -1, -1, 1, 1}; const int dj[8] = {0, 0, -1, 1, -1, 1, -1, 1}; const char dc[4] = {'U', 'D', 'L', 'R'}; bool in_grid(ll i, ll j, ll h, ll w) { return 0 <= i && i < h && 0 <= j && j < w; } #ifndef AHC_TEMPLATE_DEBUG #define AHC_TEMPLATE_DEBUG 0 #endif constexpr bool DEBUG = AHC_TEMPLATE_DEBUG; class Timer { public: using Clock = chrono::steady_clock; explicit Timer(double time_limit_seconds) : started_at_(Clock::now()), time_limit_seconds_(time_limit_seconds) {} void reset() { started_at_ = Clock::now(); } double elapsed() const { return chrono::duration(Clock::now() - started_at_).count(); } bool expired() const { return elapsed() >= time_limit_seconds_; } double time_limit() const { return time_limit_seconds_; } double remaining() const { return max(0.0, time_limit_seconds_ - elapsed()); } double progress() const { if (time_limit_seconds_ <= 0.0) return 1.0; return min(1.0, max(0.0, elapsed() / time_limit_seconds_)); } private: Clock::time_point started_at_; double time_limit_seconds_; }; // SplitMix64、高速かつシードを固定するとローカル実験を再現できる class Random { public: explicit Random(u64 seed = 0x123456789abcdef0ULL) : state_(seed) {} u64 next_u64() { u64 value = (state_ += 0x9e3779b97f4a7c15ULL); value = (value ^ (value >> 30)) * 0xbf58476d1ce4e5b9ULL; value = (value ^ (value >> 27)) * 0x94d049bb133111ebULL; return value ^ (value >> 31); } // [0, upper) int uniform_int(int upper) { assert(upper > 0); const u128 product = u128(next_u64()) * upper; return int(product >> 64); } // [lower, upper) int uniform_int(int lower, int upper) { assert(lower < upper); return lower + uniform_int(upper - lower); } // [0, 1) double uniform01() { return double(next_u64() >> 11) / 9007199254740992.0; } private: u64 state_; }; constexpr double TIME_LIMIT = 1.92; Timer timer(TIME_LIMIT); Random rng(0x123456789abcdef0ULL); struct Forest { int n; int k; vector d; vector original_d; vector danger; int total_start_weight = 0; vector cut_cells; vector> cut_lines; Forest(vector dryness, int budget) : n(int(dryness.size())), k(budget), d(move(dryness)), original_d(d), danger(n * n, 0) { rep(i, n) rep(j, n) { const int id = int(i * n + j); const int value = d[i][j] - '0'; danger[id] = value * value; rep(direction, 4) { const int ni = int(i) + di[direction]; const int nj = int(j) + dj[direction]; if (!in_grid(ni, nj, n, n)) continue; const int neighbor = d[ni][nj] - '0'; danger[id] += neighbor * neighbor; } if (value > 0) total_start_weight += danger[id]; } } bool has_tree(int id) const { return 0 <= id && id < n * n && d[id / n][id % n] != '0'; } // 抽選対象は未伐採の木のみ。重みは自身と4近傍の乾燥度の二乗和。 int choose_start(Random& random) const { if (int(cut_cells.size()) >= k || total_start_weight == 0) return -1; int ticket = random.uniform_int(total_start_weight); int id = 0; for (const string& row : d) { for (const char value : row) { if (value != '0') { if (ticket < danger[id]) return id; ticket -= danger[id]; } ++id; } } assert(false); return -1; } bool cut_tree(int id) { if (int(cut_cells.size()) >= k || !has_tree(id)) return false; const int i = id / n; const int j = id % n; const int value = d[i][j] - '0'; const int removed = value * value; total_start_weight -= danger[id]; d[i][j] = '0'; danger[id] -= removed; rep(direction, 4) { const int ni = i + di[direction]; const int nj = j + dj[direction]; if (!in_grid(ni, nj, n, n)) continue; const int neighbor_id = ni * n + nj; danger[neighbor_id] -= removed; if (has_tree(neighbor_id)) total_start_weight -= removed; } cut_cells.push_back(id); return true; } void apply_cut_line(vector line) { assert(!line.empty()); for (const int id : line) { const bool cut = cut_tree(id); assert(cut); (void)cut; } cut_lines.push_back(move(line)); } // 指定列の木を元に戻し、危険度・抽選重みを差分更新する。 void remove_cut_line(int index) { assert(0 <= index && index < int(cut_lines.size())); const vector removed_line = move(cut_lines[index]); cut_lines.erase(cut_lines.begin() + index); for (const int id : removed_line) { const int i = id / n; const int j = id % n; assert(d[i][j] == '0' && original_d[i][j] != '0'); const int value = original_d[i][j] - '0'; const int restored = value * value; d[i][j] = original_d[i][j]; danger[id] += restored; total_start_weight += danger[id]; rep(direction, 4) { const int ni = i + di[direction]; const int nj = j + dj[direction]; if (!in_grid(ni, nj, n, n)) continue; const int neighbor_id = ni * n + nj; danger[neighbor_id] += restored; if (has_tree(neighbor_id)) total_start_weight += restored; } } cut_cells.clear(); for (const auto& line : cut_lines) { cut_cells.insert(cut_cells.end(), line.begin(), line.end()); } } }; constexpr array CUT_LENGTHS = {8, 16, 24, 32}; constexpr int DRY_SIMULATION_COUNT = 3; constexpr double INITIAL_TEMPERATURE_RATIO = 0.01; constexpr double FINAL_TEMPERATURE_RATIO = 0.00001; // 方向は左上、右上、左下、右下の4種類。長さには開始地点も含める。 vector make_cut_line(const Forest& forest, int start, int direction, int length) { assert(0 <= direction && direction < 4); vector line; const int remaining = forest.k - int(forest.cut_cells.size()); const int limit = min(length, remaining); if (limit <= 0) return line; line.reserve(limit); int current = start; while (int(line.size()) < limit && forest.has_tree(current)) { line.push_back(current); if (int(line.size()) == limit) break; const int ni = current / forest.n + di[direction + 4]; const int nj = current % forest.n + dj[direction + 4]; if (!in_grid(ni, nj, forest.n, forest.n)) break; current = ni * forest.n + nj; } return line; } vector grow_cut_line(Forest& trial, int start, int direction, int length) { vector line = make_cut_line(trial, start, direction, length); for (const int id : line) { const bool cut = trial.cut_tree(id); assert(cut); (void)cut; } return line; } // floor(random * 10 / 2^64) < d と同値な閾値。 // 木がないマスでも乱数を1回消費し、従来の乱数列を保つ。 constexpr array DRY_THRESHOLDS = [] { array thresholds{}; for (int value = 1; value <= 9; ++value) { thresholds[value] = u64(((u128(value) << 64) + 9) / 10); } return thresholds; }(); vector sample_dryness(const Forest& forest, Random& random) { vector dry(forest.n * forest.n); int id = 0; for (const string& row : forest.d) { for (const char value : row) { dry[id++] = random.next_u64() < DRY_THRESHOLDS[value - '0']; } } return dry; } // 外周を0で囲み、4近傍の境界判定と除算を省く。バッファは評価間で再利用。 struct SimulationWorkspace { int n = 0; int stride = 0; vector dry; vector queue; void resize(int side) { if (side == n) return; n = side; stride = n + 2; dry.assign(stride * stride, 0); queue.resize(n * n); } void sample(const Forest& forest, Random& random) { int offset = stride + 1; for (const string& row : forest.d) { for (const char value : row) { dry[offset++] = random.next_u64() < DRY_THRESHOLDS[value - '0']; } offset += 2; } } ll square_sum() { ll score = 0; unsigned char* active = dry.data(); int* pending = queue.data(); const int end = (n + 1) * stride; for (int start = stride + 1; start < end; ++start) { if (!active[start]) continue; int head = 0; int tail = 1; pending[0] = start; active[start] = 0; while (head < tail) { const int id = pending[head++]; // 固定の4方向は展開して、反復ごとのループ制御を省く。 const int up = id - stride; if (active[up]) { active[up] = 0; pending[tail++] = up; } const int down = id + stride; if (active[down]) { active[down] = 0; pending[tail++] = down; } const int left = id - 1; if (active[left]) { active[left] = 0; pending[tail++] = left; } const int right = id + 1; if (active[right]) { active[right] = 0; pending[tail++] = right; } } score += ll(tail) * tail; } return score; } }; SimulationWorkspace& simulation_workspace(int n) { static thread_local SimulationWorkspace workspace; workspace.resize(n); return workspace; } // 乾燥した木の4近傍連結成分について、サイズの二乗和を返す。 ll component_square_sum(int n, const vector& dry) { auto& workspace = simulation_workspace(n); for (int i = 0; i < n; ++i) { copy_n(dry.data() + i * n, n, workspace.dry.data() + (i + 1) * workspace.stride + 1); } return workspace.square_sum(); } ll evaluate_forest(const Forest& forest, Random& random, const Timer* deadline = nullptr, const vector* additional_cuts = nullptr) { auto& workspace = simulation_workspace(forest.n); ll score = 0; rep(sample, DRY_SIMULATION_COUNT) { if (deadline != nullptr && deadline->expired()) return INF; workspace.sample(forest, random); if (additional_cuts != nullptr) { for (const int id : *additional_cuts) { const int padded_id = (id / forest.n + 1) * workspace.stride + id % forest.n + 1; workspace.dry[padded_id] = 0; } } score += workspace.square_sum(); } // 全候補で共通の係数1/(DRY_SIMULATION_COUNT*N*N)は比較時に省略。 return score; } struct ScoredLine { vector cells; ll score = INF; }; ScoredLine choose_cut_line(const Forest& forest, int start, Random& random, const Timer* deadline = nullptr, int fixed_length = 0) { ScoredLine best; rep(direction, 4) { const int length_count = fixed_length > 0 ? 1 : int(CUT_LENGTHS.size()); rep(length_index, length_count) { const int length = fixed_length > 0 ? fixed_length : CUT_LENGTHS[length_index]; if (deadline != nullptr && deadline->expired()) return best; vector line = make_cut_line(forest, start, int(direction), length); const ll score = evaluate_forest(forest, random, deadline, &line); // 同点では先に調べた候補を保持する。 if (score < best.score) { best.score = score; best.cells = move(line); } } } return best; } void solve(Forest& forest, Random& random) { while (true) { const int start = forest.choose_start(random); if (start == -1) break; ScoredLine best = choose_cut_line(forest, start, random); forest.apply_cut_line(move(best.cells)); } } array sample_split_lengths(int length, Random& random) { assert(length >= 10); const int center = length / 2; const int radius = length / 10 + 1; return {random.uniform_int(center - radius, center + radius + 1), random.uniform_int(center - radius, center + radius + 1)}; } // 開始地点を危険度の重み付き抽選で選び、最良の列を設置する。 // fixed_length == 0 のときは初期解と同じ4種類の長さを試す。 ll place_best_line(Forest& forest, Random& random, const Timer* deadline = nullptr, int fixed_length = 0) { if (deadline != nullptr && deadline->expired()) return INF; const int start = forest.choose_start(random); if (start == -1) return INF; ScoredLine best = choose_cut_line(forest, start, random, deadline, fixed_length); if (deadline != nullptr && deadline->expired()) return INF; if (best.cells.empty()) return INF; forest.apply_cut_line(move(best.cells)); return best.score; } struct ExponentialTemperature { double initial_temperature; double log_cooling_ratio; explicit ExponentialTemperature(ll initial_score) : initial_temperature(double(initial_score) * INITIAL_TEMPERATURE_RATIO), log_cooling_ratio(log(FINAL_TEMPERATURE_RATIO / INITIAL_TEMPERATURE_RATIO)) {} double at(double progress) const { if (initial_temperature <= 0.0) return 0.0; progress = clamp(progress, 0.0, 1.0); return initial_temperature * exp(log_cooling_ratio * progress); } }; bool accept_annealing_move(ll current_score, ll candidate_score, double temperature, Random& random) { if (candidate_score == INF) return false; if (candidate_score <= current_score) return true; if (temperature <= 0.0) return false; const double delta = double(candidate_score - current_score); return random.uniform01() < exp(-delta / temperature); } // 既存の伐採位置から伸ばせる方向が残っているかを調べる。 bool has_cut_extension(const Forest& forest) { for (const int anchor : forest.cut_cells) { const int i = anchor / forest.n; const int j = anchor % forest.n; for (int direction = 4; direction < 8; ++direction) { const int ni = i + di[direction]; const int nj = j + dj[direction]; if (in_grid(ni, nj, forest.n, forest.n) && forest.has_tree(ni * forest.n + nj)) { return true; } } } return false; } int fill_remaining_cuts(Forest& forest, Random& random, const Timer* deadline = nullptr) { const int initial_count = int(forest.cut_cells.size()); while (int(forest.cut_cells.size()) < forest.k && !forest.cut_cells.empty()) { if (deadline != nullptr && deadline->expired()) break; const int anchor = forest.cut_cells[random.uniform_int(int(forest.cut_cells.size()))]; const int direction = random.uniform_int(4); const int ni = anchor / forest.n + di[direction + 4]; const int nj = anchor % forest.n + dj[direction + 4]; if (!in_grid(ni, nj, forest.n, forest.n) || !forest.has_tree(ni * forest.n + nj)) { if (!has_cut_extension(forest)) break; continue; } const int remaining = forest.k - int(forest.cut_cells.size()); forest.apply_cut_line(make_cut_line(forest, ni * forest.n + nj, direction, remaining)); } return int(forest.cut_cells.size()) - initial_count; } bool try_relocate_line(Forest& forest, int index, ll& current_score, Random& random, const Timer* deadline = nullptr, double temperature = 0.0) { if (deadline != nullptr && deadline->expired()) return false; const int removed_length = int(forest.cut_lines[index].size()); Forest base = forest; base.remove_cut_line(index); Forest candidate = base; ll candidate_score = place_best_line(candidate, random, deadline); if (removed_length >= 10) { if (deadline != nullptr && deadline->expired()) return false; const auto lengths = sample_split_lengths(removed_length, random); Forest split_candidate = base; const ll first_score = place_best_line(split_candidate, random, deadline, lengths[0]); if (first_score == INF) return false; const ll split_score = place_best_line(split_candidate, random, deadline, lengths[1]); if (split_score < candidate_score) { candidate_score = split_score; candidate = move(split_candidate); } } if (deadline != nullptr && deadline->expired()) return false; if (fill_remaining_cuts(candidate, random, deadline) > 0) { candidate_score = evaluate_forest(candidate, random, deadline); } if (deadline != nullptr && deadline->expired()) return false; if (!accept_annealing_move(current_score, candidate_score, temperature, random)) return false; forest = move(candidate); current_score = candidate_score; return true; } ll simulated_annealing(Forest& forest, Random& random, const Timer& deadline) { if (forest.cut_lines.empty() || deadline.expired()) return INF; ll current_score = evaluate_forest(forest, random, &deadline); if (current_score == INF) return INF; const ll initial_score = current_score; ll best_score = current_score; Forest best_forest = forest; const double started_at = deadline.elapsed(); const double duration = deadline.time_limit() - started_at; if (duration <= 0.0) return best_score; const double inverse_duration = 1.0 / duration; const ExponentialTemperature temperature_schedule(initial_score); u64 iterations = 0; u64 accepted = 0; while (true) { const double now = deadline.elapsed(); if (now >= deadline.time_limit()) break; ++iterations; const double progress = (now - started_at) * inverse_duration; const double temperature = temperature_schedule.at(progress); const int index = random.uniform_int(int(forest.cut_lines.size())); if (try_relocate_line(forest, index, current_score, random, &deadline, temperature)) { ++accepted; if (current_score < best_score) { best_score = current_score; best_forest = forest; } } } forest = move(best_forest); if constexpr (DEBUG) { cerr << "[sa] iterations=" << iterations << " accepted=" << accepted << " elapsed=" << deadline.elapsed() << '\n'; } return best_score; } #ifndef AHC_TEMPLATE_NO_MAIN int main() { ios::sync_with_stdio(false); cin.tie(nullptr); int n, k; if (!(cin >> n >> k)) return 0; vector d(n); for (string& row : d) cin >> row; Forest forest(move(d), k); solve(forest, rng); simulated_annealing(forest, rng, timer); cout << forest.cut_cells.size() << '\n'; for (const int id : forest.cut_cells) { cout << id / n << ' ' << id % n << '\n'; } return 0; } #endif