#include #pragma GCC target("avx2") #pragma GCC optimize("O3") #pragma GCC optimize("unroll-loops") #include using namespace std; // MT19937 のパラメータ・出力列を保ち、状態の格納幅だけを 32 bit にする。 using SamplingRng = mersenne_twister_engine; using lint = long long int; using P = pair; using PL = pair; #define FOR(i, begin, end) for(int i=(begin),i##_end_=(end);i=i##_begin_;i--) #define REP(i, n) FOR(i,0,n) #define IREP(i, n) IFOR(i,0,n) #define ALL(a) (a).begin(),(a).end() constexpr int MOD = 1000000007; vector RH_B = {1532834020, 1388622299}; vector RH_M = {2147482409, 2147478017}; constexpr int INF = 2147483647; void yes(bool expr) {cout << (expr ? "Yes" : "No") << "\n";} templatevoid chmax(T &a, const T &b) { if (avoid chmin(T &a, const T &b) { if (b> neighbors; explicit GridGraph(int side) : n(side), cells(n * n), neighbors(cells) { constexpr int dr[4] = {-1, 1, 0, 0}; constexpr int dc[4] = {0, 0, -1, 1}; for (int r = 0; r < n; ++r) { for (int c = 0; c < n; ++c) { auto &adj = neighbors[r * n + c]; adj.fill(-1); for (int direction = 0; direction < 4; ++direction) { int nr = r + dr[direction]; int nc = c + dc[direction]; if (0 <= nr && nr < n && 0 <= nc && nc < n) { adj[direction] = nr * n + nc; } } } } } }; // N=50 の制約では頂点番号・成分サイズが 16 bit に収まる。 // 外周を番兵で囲み、order=0 を木なし、-1 を未訪問の乾燥木とする。 struct ArticulationWorkspace { static constexpr int MAX_SIDE = 50; static constexpr int CAPACITY = (MAX_SIDE + 2) * (MAX_SIDE + 2); array positions; array order{}; explicit ArticulationWorkspace(int side) : cells(side * side) { const int stride = side + 2; offsets = {-stride, stride, -1, 1}; for (int r = 0; r < side; ++r) { for (int c = 0; c < side; ++c) { const int flat = r * side + c; const int padded = (r + 1) * stride + c + 1; positions[flat] = static_cast(padded); flat_index[padded] = static_cast(flat); } } } void accumulate(vector &total_decrease) { timer = 0; for (int flat = 0; flat < cells; ++flat) { const int root = positions[flat]; if (order[root] != -1) continue; component_count = 0; dfs(root, -1); const int size = subtree_size[root]; const long long square = 1LL * size * size; for (int i = 0; i < component_count; ++i) { const int v = component[i]; if (!articulation[v]) continue; const int rest = size - 1 - separated_size[v]; total_decrease[flat_index[v]] += square - separated_squares[v] - 1LL * rest * rest; } } } private: int cells; array offsets; array flat_index, low, subtree_size, separated_size, component; array separated_squares; array articulation; int timer = 0, component_count = 0; void dfs(int v, int parent) { const int own_order = ++timer; order[v] = static_cast(own_order); int lowest = own_order, size = 1, separated = 0, parts = 0, squares = 0; for (int offset : offsets) { const int next = v + offset; const int next_order = order[next]; if (next_order == 0) continue; if (next_order == -1) { dfs(next, v); const int child_size = subtree_size[next]; const int child_low = low[next]; size += child_size; lowest = min(lowest, child_low); if (child_low >= own_order) { separated += child_size; squares += child_size * child_size; ++parts; } } else if (next != parent) { lowest = min(lowest, next_order); } } low[v] = static_cast(lowest); subtree_size[v] = static_cast(size); separated_size[v] = static_cast(separated); separated_squares[v] = squares; articulation[v] = parent == -1 ? parts >= 2 : parts >= 1; component[component_count++] = static_cast(v); } }; // 提出時のサンプリング回数。実行時は --samples X でも変更できる。 constexpr int DEFAULT_SAMPLE_COUNT = 1200; constexpr int SCORE_SAMPLE_COUNT = 100; constexpr int CUT_BATCH_SIZE = 3; constexpr int DRYNESS_CORRECTION_START_PERCENT = 70; constexpr int DRYNESS_CORRECTION_END_PERCENT = 100; long long component_square_sum(const GridGraph &graph, const vector &dry) { vector visited(graph.cells, false); vector component; component.reserve(graph.cells); long long total = 0; for (int root = 0; root < graph.cells; ++root) { if (!dry[root] || visited[root]) continue; component.clear(); component.push_back(root); visited[root] = true; for (size_t head = 0; head < component.size(); ++head) { for (int next : graph.neighbors[component[head]]) { if (next == -1 || !dry[next] || visited[next]) continue; visited[next] = true; component.push_back(next); } } const long long size = static_cast(component.size()); total += size * size; } return total; } vector select_best_cuts(const vector &decrease, const vector &available, int max_count, int side) { if (max_count <= 0) return {}; if (max_count == 1) { int best = -1; for (int v = 0; v < static_cast(available.size()); ++v) { if (available[v] && (best == -1 || decrease[v] > decrease[best])) best = v; } if (best == -1) return {}; return {best}; } vector candidates; candidates.reserve(available.size()); for (int v = 0; v < static_cast(available.size()); ++v) { if (available[v]) candidates.push_back(v); } const int count = min(max_count, static_cast(candidates.size())); const auto lower_priority = [&](int a, int b) { if (decrease[a] != decrease[b]) return decrease[a] < decrease[b]; return a > b; }; make_heap(candidates.begin(), candidates.end(), lower_priority); vector blocked(available.size(), false); vector selected; selected.reserve(count); while (!candidates.empty() && static_cast(selected.size()) < count) { pop_heap(candidates.begin(), candidates.end(), lower_priority); const int v = candidates.back(); candidates.pop_back(); if (blocked[v]) continue; selected.push_back(v); // この step で選んだ木からマンハッタン距離 3 以内を除外する。 const int row = v / side, col = v % side; for (int r = max(0, row - 12); r <= min(side - 1, row + 12); ++r) { const int radius = 12 - abs(r - row); const int first = r * side + max(0, col - radius); const int last = r * side + min(side - 1, col + radius) + 1; fill(blocked.begin() + first, blocked.begin() + last, true); } } return selected; } int main(int argc, char **argv) { ios::sync_with_stdio(false); cin.tie(nullptr); int sample_count = DEFAULT_SAMPLE_COUNT; if (argc != 1) { if (argc != 3 || string(argv[1]) != "--samples") { cerr << "Usage: " << argv[0] << " [--samples X]\n"; return 1; } const char *begin = argv[2]; const char *end = begin + strlen(begin); auto parsed = from_chars(begin, end, sample_count); if (parsed.ec != errc{} || parsed.ptr != end || sample_count <= 0) { cerr << "X must be a positive integer.\n"; return 1; } } int n, k; if (!(cin >> n >> k)) return 0; vector board(n); for (auto &row : board) cin >> row; GridGraph graph(n); vector dryness_by_cell(graph.cells); vector available(graph.cells); for (int v = 0; v < graph.cells; ++v) { dryness_by_cell[v] = static_cast(board[v / n][v % n] - '0'); available[v] = dryness_by_cell[v] != 0; } // 乱数生成器は一度だけ初期化し、各 step で新しい乱数を取得する。 SamplingRng rng(0); uniform_int_distribution random_digit(0, 9); uniform_int_distribution random_digit_100(0, 99); ArticulationWorkspace sample(n); const int batch_count = (k + CUT_BATCH_SIZE - 1) / CUT_BATCH_SIZE; const long long correction_steps = max(1, batch_count - 1); const long long threshold_denominator = 10 * correction_steps; vector answer; vector total_decrease(graph.cells); for (int step = 0; step < batch_count; ++step) { // 補正係数を初回から最終回まで線形に変化させる。 // 整数の比で閾値を表し、境界での浮動小数点誤差を避ける。 const long long correction_numerator = DRYNESS_CORRECTION_START_PERCENT * correction_steps + 1LL * (DRYNESS_CORRECTION_END_PERCENT - DRYNESS_CORRECTION_START_PERCENT) * step; array dry_cutoff; for (int d = 0; d <= 9; ++d) { // R * 分母 < d * 分子 と、R < ceil(d * 分子 / 分母) は同値。 dry_cutoff[d] = static_cast((d * correction_numerator + threshold_denominator - 1) / threshold_denominator); } fill(total_decrease.begin(), total_decrease.end(), 0); for (int s = 0; s < sample_count; ++s) { // 乾燥状態を毎回生成し直し、伐採済みの木を除外する。 for (int v = 0; v < graph.cells; ++v) { const int random_value = random_digit_100(rng); sample.order[sample.positions[v]] = available[v] && random_value < dry_cutoff[dryness_by_cell[v]] ? -1 : 0; } sample.accumulate(total_decrease); } // 評価値順に、同じ step 内の距離制約を満たす M 本を選ぶ。 // 同点は座標順。 const int remaining = k - static_cast(answer.size()); const auto selected = select_best_cuts(total_decrease, available, min(CUT_BATCH_SIZE, remaining), n); if (selected.empty()) break; for (int v : selected) { answer.push_back(v); available[v] = false; } } cout << answer.size() << '\n'; for (int v : answer) cout << v / n << ' ' << v % n << '\n'; // 伐採の選択に使ったサンプルとは別に、100 個を新しく生成する。 // 各サンプルを伐採前・伐採後で共用して E0 と E を推定する。 long long before_sum = 0, after_sum = 0; for (int s = 0; s < SCORE_SAMPLE_COUNT; ++s) { vector dry(graph.cells); for (int v = 0; v < graph.cells; ++v) { dry[v] = random_digit(rng) < dryness_by_cell[v]; } before_sum += component_square_sum(graph, dry); for (int v : answer) dry[v] = false; after_sum += component_square_sum(graph, dry); } const long double denominator = 1.0L * SCORE_SAMPLE_COUNT * graph.cells; const long double e0 = before_sum / denominator; const long double e = after_sum / denominator; const long double score = round(1.0e6L * (1.0L + log2(e0 / e))); cerr << fixed << setprecision(0) << score << '\n'; return 0; }