#include #include #include #include #include #include #include #include #include #include #include #include using namespace std; using Connection = pair; constexpr int SAFE_ZONE_CRITERIA = 4; constexpr int CAND_SIZES_LIM = 30000; constexpr double TIME_LIMIT = 1.8; constexpr double ROUTE_TIME_LIMIT = 1.0; constexpr double MINCUT_TIME_LIMIT = 0.18; constexpr double MINCUT_PROB = 0.5; constexpr int REFINE_SAMPLES = 512; constexpr double START_TEMP = 10000; constexpr double END_TEMP = 100; constexpr double TWO_ADD_ONE_REMOVE_PROB = 0.1; constexpr double TWO_REMOVE_ONE_ADD_PROB = 0.1; constexpr double SHUFFLE_PROB = 0.1; constexpr double TWO_ADD_TWO_REMOVE_PROB = 0.1; int N, K, N2; vector> board; vector danger_zones; int to_1d(int i, int j) { return i * N2 + j; } pair to_2d(int v) { return {v / N2, v % N2}; } int connection_dist(int v1, int v2) { auto [x1, y1] = to_2d(v1); auto [x2, y2] = to_2d(v2); return min(abs(x1 - x2), abs(y1 - y2)); } int connection_dist(Connection edge) { return connection_dist(edge.first, edge.second); } vector create_route(int v1, int v2) { auto [x1, y1] = to_2d(v1); auto [x2, y2] = to_2d(v2); if (x1 > x2) { swap(x1, x2); swap(y1, y2); } vector route; int x = x1, y = y1; while (x != x2 || y != y2) { if (x < x2) ++x; if (y < y2) ++y; else if (y > y2) --y; if (x != x2 || y != y2) route.push_back(to_1d(x, y)); } return route; } struct Dinic { struct Edge { int to, next, capacity; }; vector head, level, cursor, queue; vector edges; explicit Dinic(int n) : head(n, -1), level(n), cursor(n), queue(n) { edges.reserve(n * 6); } void add_edge(int from, int to, int capacity) { edges.push_back({to, head[from], capacity}); head[from] = static_cast(edges.size()) - 1; edges.push_back({from, head[to], 0}); head[to] = static_cast(edges.size()) - 1; } bool bfs(int source, int sink) { fill(level.begin(), level.end(), -1); int begin = 0, end = 0; queue[end++] = source; level[source] = 0; while (begin < end) { int v = queue[begin++]; for (int e = head[v]; e >= 0; e = edges[e].next) { const auto& edge = edges[e]; if (edge.capacity > 0 && level[edge.to] < 0) { level[edge.to] = level[v] + 1; queue[end++] = edge.to; } } } return level[sink] >= 0; } int dfs(int v, int sink, int flow) { if (v == sink) return flow; for (int& e = cursor[v]; e >= 0; e = edges[e].next) { auto& edge = edges[e]; if (edge.capacity == 0 || level[edge.to] != level[v] + 1) continue; int sent = dfs(edge.to, sink, min(flow, edge.capacity)); if (sent) { edge.capacity -= sent; edges[e ^ 1].capacity += sent; return sent; } } return 0; } int max_flow(int source, int sink, int limit) { int total = 0; while (total < limit && bfs(source, sink)) { cursor = head; while (total < limit) { int sent = dfs(source, sink, limit - total); if (!sent) break; total += sent; } } return total; } }; // 4近傍の危険領域を頂点カットで分割する。各木を in -> out に分け、 // この辺に伐採コスト、隣接辺に十分大きい容量を置く。 // 射影の両端を端子領域にして、1本だけを孤立させるカットへの偏りを抑える。 vector> generate_mincuts(uint32_t seed) { const auto start = chrono::steady_clock::now(); vector> components, result; vector label(N2 * N2, -1), local(N2 * N2, -1); for (int v : danger_zones) if (label[v] < 0) { int id = static_cast(components.size()); components.push_back({v}); label[v] = id; for (int j = 0; j < static_cast(components.back().size()); ++j) { int u = components.back()[j]; for (int d : {-N2, N2, -1, 1}) { int w = u + d; if (board[w / N2][w % N2] > SAFE_ZONE_CRITERIA && label[w] < 0) { label[w] = id; components.back().push_back(w); } } } } sort(components.begin(), components.end(), [](const auto& a, const auto& b) { return a.size() > b.size(); }); while (!components.empty() && components.back().size() < 40) components.pop_back(); if (components.empty()) return result; mt19937 rng(seed ^ 0x9e3779b9U); int mass = 0; for (const auto& component : components) mass += static_cast(component.size() * component.size()); constexpr array, 8> directions = {{{1, 0}, {0, 1}, {1, 1}, {1, -1}, {2, 1}, {2, -1}, {1, 2}, {1, -2}}}; constexpr int INF = 100000000; for (int trial = 0; trial < 256; ++trial) { if (chrono::duration(chrono::steady_clock::now() - start).count() >= MINCUT_TIME_LIMIT) break; int pick = static_cast(rng() % mass), id = 0; while (pick >= static_cast(components[id].size() * components[id].size())) { pick -= static_cast(components[id].size() * components[id].size()); ++id; } const auto& vertices = components[id]; int n = static_cast(vertices.size()); auto [dx, dy] = directions[trial % directions.size()]; vector projection(n), sorted_projection; for (int j = 0; j < n; ++j) { int v = vertices[j]; local[v] = j; projection[j] = dx * (v / N2) + dy * (v % N2); } sorted_projection = projection; sort(sorted_projection.begin(), sorted_projection.end()); int center = 20 + static_cast(rng() % 61); int half_width = 10 + static_cast(rng() % 21); int low = sorted_projection[n * max(5, center - half_width) / 100]; int high = sorted_projection[n * min(95, center + half_width) / 100]; if (low >= high) continue; int source = 2 * n, sink = source + 1; Dinic graph(2 * n + 2); for (int j = 0; j < n; ++j) { bool terminal = projection[j] <= low || projection[j] >= high; // 小さな容量の揺らぎで、ほぼ同じ伐採本数の異なるカットを得る。 graph.add_edge(2 * j, 2 * j + 1, terminal ? INF : 1000 + static_cast(rng() % 300)); if (projection[j] <= low) graph.add_edge(source, 2 * j, INF); if (projection[j] >= high) graph.add_edge(2 * j + 1, sink, INF); int v = vertices[j]; for (int d : {-N2, N2, -1, 1}) { int u = v + d; // 別の連結成分には辺がなく、境界の番兵マスは常に湿っている。 if (board[u / N2][u % N2] > SAFE_ZONE_CRITERIA) { graph.add_edge(2 * j + 1, 2 * local[u], INF); } } } if (graph.max_flow(source, sink, 1300 * K + 1) >= 1300 * K + 1) continue; vector cut; for (int j = 0; j < n; ++j) { if (graph.level[2 * j] >= 0 && graph.level[2 * j + 1] < 0) cut.push_back(vertices[j]); } if (cut.empty() || static_cast(cut.size()) > K) continue; sort(cut.begin(), cut.end()); if (find(result.begin(), result.end(), cut) == result.end()) result.push_back(std::move(cut)); } return result; } // 入力と候補辺だけで決まる情報をキャッシュし、反復用メモリを再利用する。 struct SearchWorkspace { vector routes, route_offsets; vector distances, danger_idx, trees, parents; vector> adjacent_dangers; vector used; vector removal_mask, cached_removal_mask; vector removed_indices; long long cached_score = 0; bool score_cached = false; unsigned int tree_stamp = 0; void add_cut(const vector& cut) { routes.insert(routes.end(), cut.begin(), cut.end()); route_offsets.push_back(static_cast(routes.size())); distances.push_back(static_cast(cut.size())); } explicit SearchWorkspace(const vector& candidates) : danger_idx(N2 * N2, -1), parents(danger_zones.size()), used(N2 * N2, 0), removal_mask((danger_zones.size() + 63) / 64), cached_removal_mask(removal_mask.size()) { routes.reserve(candidates.size() * N2); route_offsets.reserve(candidates.size() + 1); route_offsets.push_back(0); distances.reserve(candidates.size()); trees.reserve(N2 * N2); removed_indices.reserve(N2 * N2); for (auto [v1, v2] : candidates) { auto [x1, y1] = to_2d(v1); auto [x2, y2] = to_2d(v2); if (x1 > x2) { swap(x1, x2); swap(y1, y2); } int x = x1, y = y1; while (x != x2 || y != y2) { if (x < x2) ++x; if (y < y2) ++y; else if (y > y2) --y; if ((x != x2 || y != y2) && board[x][y] != 0) { routes.push_back(to_1d(x, y)); } } route_offsets.push_back(static_cast(routes.size())); distances.push_back(connection_dist(v1, v2)); } for (int i = 0; i < static_cast(danger_zones.size()); ++i) { danger_idx[danger_zones[i]] = i; } for (int i = 0; i < static_cast(danger_zones.size()); ++i) { int v = danger_zones[i]; for (int d : {1, N2}) { int neighbor = v + d; if (neighbor < N2 * N2 && danger_idx[neighbor] >= 0) { adjacent_dangers.emplace_back(i, danger_idx[neighbor]); } } } } const vector& build(const vector& edges) { if (++tree_stamp == 0) { fill(used.begin(), used.end(), 0); tree_stamp = 1; } trees.clear(); for (auto it = edges.rbegin(); it != edges.rend(); ++it) { for (int i = route_offsets[*it]; i < route_offsets[*it + 1]; ++i) { int v = routes[i]; if (used[v] == tree_stamp) continue; used[v] = tree_stamp; trees.push_back(v); if (static_cast(trees.size()) == K) return trees; } } return trees; } int find_root(int x) { if (parents[x] < 0) return x; return parents[x] = find_root(parents[x]); } long long evaluate(const vector& selected_trees) { constexpr int REMOVED = numeric_limits::min(); fill(removal_mask.begin(), removal_mask.end(), 0); removed_indices.clear(); for (int v : selected_trees) { int idx = danger_idx[v]; if (idx < 0) continue; uint64_t bit = uint64_t{1} << (idx % 64); if ((removal_mask[idx / 64] & bit) == 0) { removal_mask[idx / 64] |= bit; removed_indices.push_back(idx); } } // 削除対象が前回と完全一致する場合だけ、同じ評価値を再利用する。 if (score_cached && removal_mask == cached_removal_mask) return cached_score; fill(parents.begin(), parents.end(), -1); for (int idx : removed_indices) parents[idx] = REMOVED; // 無向辺は一度だけ処理する。連結成分と評価値は元の処理と同じ。 for (auto [a, b] : adjacent_dangers) { if (parents[a] == REMOVED || parents[b] == REMOVED) continue; a = find_root(a); b = find_root(b); if (a == b) continue; if (parents[a] > parents[b]) swap(a, b); parents[a] += parents[b]; parents[b] = a; } long long score = 0; for (int size : parents) { if (size < 0 && size != REMOVED) score += 1LL * size * size; } cached_removal_mask = removal_mask; cached_score = score; score_cached = true; return score; } }; // 共通の乾燥サンプルで比較する。各頂点を削除したときの改善量は // DFS 木の low-link と部分木サイズから、全頂点分を線形時間で求める。 struct SampleEvaluator { static constexpr int SAMPLES = 128; int cells; vector vertices; vector> neighbors; vector> dry; vector cut; vector tin, low, subtree, component, component_sizes, order, separated; vector separated_sq, gain, restore_cost; int timer = 0; explicit SampleEvaluator(int sample_count = SAMPLES) : cells(N2 * N2), neighbors(cells), dry(sample_count, vector(cells)), cut(cells), tin(cells), low(cells), subtree(cells), component(cells), separated(cells), separated_sq(cells), gain(cells), restore_cost(cells) { uint64_t seed = 1469598103934665603ULL; for (int x = 1; x <= N; ++x) for (int y = 1; y <= N; ++y) { int v = to_1d(x, y); neighbors[v] = {v - N2, v + N2, v - 1, v + 1}; seed = (seed ^ uint64_t(board[x][y])) * 1099511628211ULL; if (board[x][y]) vertices.push_back(v); } mt19937_64 random(seed); for (auto& sample : dry) for (int v : vertices) { sample[v] = random() % 10 < uint64_t(board[v / N2][v % N2]); } order.reserve(cells); component_sizes.reserve(cells); } void dfs(int v, int parent, int id, const vector& active) { tin[v] = low[v] = ++timer; subtree[v] = 1; separated[v] = 0; separated_sq[v] = 0; component[v] = id; order.push_back(v); for (int u : neighbors[v]) { if (!active[u] || cut[u] || u == parent) continue; if (tin[u]) { low[v] = min(low[v], tin[u]); } else { dfs(u, v, id, active); subtree[v] += subtree[u]; low[v] = min(low[v], low[u]); if (low[u] >= tin[v]) { separated[v] += subtree[u]; separated_sq[v] += 1LL * subtree[u] * subtree[u]; } } } } long long analyze() { fill(gain.begin(), gain.end(), 0); fill(restore_cost.begin(), restore_cost.end(), 0); long long score = 0; for (const auto& sample : dry) { fill(tin.begin(), tin.end(), 0); component_sizes.clear(); order.clear(); timer = 0; for (int v : vertices) if (sample[v] && !cut[v] && !tin[v]) { int id = static_cast(component_sizes.size()); int begin = static_cast(order.size()); dfs(v, -1, id, sample); int size = subtree[v]; component_sizes.push_back(size); score += 1LL * size * size; for (int j = begin; j < static_cast(order.size()); ++j) { int u = order[j]; int rest = size - 1 - separated[u]; gain[u] += 1LL * size * size - separated_sq[u] - 1LL * rest * rest; } } for (int v : vertices) if (cut[v] && sample[v]) { array ids; int count = 0, sum = 1; long long square_sum = 0; for (int u : neighbors[v]) if (sample[u] && !cut[u]) { int id = component[u]; bool duplicate = false; for (int j = 0; j < count; ++j) duplicate |= ids[j] == id; if (duplicate) continue; ids[count++] = id; int size = component_sizes[id]; sum += size; square_sum += 1LL * size * size; } restore_cost[v] += 1LL * sum * sum - square_sum; } } return score; } }; // 変更した木を含む成分と、復元した木に隣接する成分だけを作り直す。 // サンプルごとの頂点削除効果を保持し、全サンプルの和を差分更新する。 struct DynamicEvaluator { struct State { vector component, next, size, gain; explicit State(int n) : component(n, -1), next(n, -1), size(n), gain(n) {} }; const SampleEvaluator& samples; vector states; vector cut; vector gain, restore_cost; vector tin, low, subtree, separated, separated_sq, order, affected, roots, inserted; int timer = 0; long long score = 0; DynamicEvaluator(const SampleEvaluator& samples, const vector& trees) : samples(samples), cut(samples.cells), gain(samples.cells), restore_cost(samples.cells), tin(samples.cells), low(samples.cells), subtree(samples.cells), separated(samples.cells), separated_sq(samples.cells) { for (int v : trees) cut[v] = 1; states.reserve(samples.dry.size()); order.reserve(samples.cells); affected.reserve(samples.cells); for (int s = 0; s < static_cast(samples.dry.size()); ++s) { states.emplace_back(samples.cells); fill(tin.begin(), tin.end(), 0); timer = 0; for (int v : samples.vertices) if (samples.dry[s][v] && !cut[v] && !tin[v]) { build_component(s, v); } } } void dfs(int s, int v, int parent, int root) { auto& state = states[s]; tin[v] = low[v] = ++timer; subtree[v] = 1; separated[v] = separated_sq[v] = 0; state.component[v] = root; order.push_back(v); for (int u : samples.neighbors[v]) { if (!samples.dry[s][u] || cut[u] || u == parent) continue; if (tin[u]) { low[v] = min(low[v], tin[u]); } else { dfs(s, u, v, root); subtree[v] += subtree[u]; low[v] = min(low[v], low[u]); if (low[u] >= tin[v]) { separated[v] += subtree[u]; separated_sq[v] += subtree[u] * subtree[u]; } } } } void build_component(int s, int root) { auto& state = states[s]; order.clear(); dfs(s, root, -1, root); int size = static_cast(order.size()); state.size[root] = size; score += 1LL * size * size; for (int i = 0; i < size; ++i) { int v = order[i]; state.next[v] = i + 1 < size ? order[i + 1] : -1; int rest = size - 1 - separated[v]; state.gain[v] = size * size - separated_sq[v] - rest * rest; gain[v] += state.gain[v]; } } void invalidate(int s, int root) { auto& state = states[s]; score -= 1LL * state.size[root] * state.size[root]; state.size[root] = 0; for (int v = root; v >= 0; v = state.next[v]) { gain[v] -= state.gain[v]; state.gain[v] = 0; state.component[v] = -1; tin[v] = 0; affected.push_back(v); } } // changes は相異なる木と、その変更後の伐採状態の組。 void apply(const vector>& changes) { for (auto [v, removed] : changes) cut[v] = removed; for (int s = 0; s < static_cast(states.size()); ++s) { auto& state = states[s]; roots.clear(); inserted.clear(); for (auto [v, removed] : changes) { if (!samples.dry[s][v]) continue; if (state.component[v] >= 0) roots.push_back(state.component[v]); else if (!removed) inserted.push_back(v); if (!removed) for (int u : samples.neighbors[v]) { if (state.component[u] >= 0) roots.push_back(state.component[u]); } } if (roots.empty() && inserted.empty()) continue; sort(roots.begin(), roots.end()); roots.erase(unique(roots.begin(), roots.end()), roots.end()); affected.clear(); for (int root : roots) invalidate(s, root); for (int v : inserted) { tin[v] = 0; affected.push_back(v); } timer = 0; for (int v : affected) if (!cut[v] && !tin[v]) build_component(s, v); } } void restoration_costs(const vector& trees) { fill(restore_cost.begin(), restore_cost.end(), 0); for (int s = 0; s < static_cast(states.size()); ++s) { const auto& state = states[s]; for (int v : trees) if (samples.dry[s][v]) { array ids; int count = 0, sum = 1; long long square_sum = 0; for (int u : samples.neighbors[v]) { int root = state.component[u]; if (root < 0) continue; bool duplicate = false; for (int j = 0; j < count; ++j) duplicate |= ids[j] == root; if (duplicate) continue; ids[count++] = root; int size = state.size[root]; sum += size; square_sum += 1LL * size * size; } restore_cost[v] += 1LL * sum * sum - square_sum; } } } }; void refine_trees(vector& trees, chrono::steady_clock::time_point start) { SampleEvaluator samples(REFINE_SAMPLES); DynamicEvaluator evaluator(samples, trees); vector tried(samples.cells); auto expired = [&]() { return chrono::duration(chrono::steady_clock::now() - start).count() >= TIME_LIMIT; }; int iterations = 0, improvements = 0; bool costs_dirty = true; while (!expired()) { if (static_cast(trees.size()) < K) { int add = -1; for (int v : samples.vertices) if (!evaluator.cut[v] && (add < 0 || evaluator.gain[v] > evaluator.gain[add])) add = v; if (add < 0) break; evaluator.apply({{add, true}}); trees.push_back(add); costs_dirty = true; } else { // 棄却した交換では状態が元に戻るため、復元コストも再利用できる。 if (costs_dirty) evaluator.restoration_costs(trees); costs_dirty = false; int remove = -1; for (int v : trees) if (!tried[v] && (remove < 0 || evaluator.restore_cost[v] < evaluator.restore_cost[remove])) remove = v; if (remove < 0) break; long long previous_score = evaluator.score; evaluator.apply({{remove, false}}); int add = remove; for (int v : samples.vertices) if (!evaluator.cut[v] && evaluator.gain[v] > evaluator.gain[add]) add = v; if (evaluator.score - evaluator.gain[add] < previous_score) { evaluator.apply({{add, true}}); *find(trees.begin(), trees.end(), remove) = add; fill(tried.begin(), tried.end(), 0); ++improvements; costs_dirty = true; } else { evaluator.apply({{remove, true}}); tried[remove] = 1; } } ++iterations; } cerr << "dynamic_iterations=" << iterations << " improvements=" << improvements << " sample_score=" << evaluator.score << '\n'; } int main() { const auto start_time = chrono::steady_clock::now(); ios::sync_with_stdio(false); cin.tie(nullptr); if (!(cin >> N >> K)) return 0; N2 = N + 2; board.assign(N2, vector(N2, 0)); for (int i = 0; i < N; ++i) { string row; cin >> row; for (int j = 0; j < N; ++j) board[i + 1][j + 1] = row[j] - '0'; } uint32_t input_seed = 2166136261U; for (const auto& row : board) for (int d : row) { input_seed = (input_seed ^ static_cast(d)) * 16777619U; } mt19937 rng(input_seed); auto randint = [&](int lo, int hi) { return uniform_int_distribution(lo, hi)(rng); }; auto random_real = [&]() { return uniform_real_distribution(0.0, 1.0)(rng); }; vector safe_zones; for (int i = 0; i < N2; ++i) { for (int j = 0; j < N2; ++j) { if (board[i][j] <= SAFE_ZONE_CRITERIA) safe_zones.push_back(to_1d(i, j)); else danger_zones.push_back(to_1d(i, j)); } } unordered_set candidate_set; candidate_set.reserve(CAND_SIZES_LIM + 1); vector edge_candidates; edge_candidates.reserve(CAND_SIZES_LIM + 1); // 候補総数が上限より少ない入力でも終了するよう、試行回数を制限する。 for (int attempts = 0; candidate_set.size() <= CAND_SIZES_LIM && attempts < 1000000; ++attempts) { int v1 = safe_zones[randint(0, static_cast(safe_zones.size()) - 1)]; int v2 = safe_zones[randint(0, static_cast(safe_zones.size()) - 1)]; int dist = connection_dist(v1, v2); if (2 <= dist && dist <= 50) { uint64_t key = (uint64_t(static_cast(v1)) << 32) | static_cast(v2); if (candidate_set.insert(key).second) edge_candidates.emplace_back(v1, v2); } } // 旧版の set と同じ候補順に戻し、乱数に対する選択結果を維持する。 sort(edge_candidates.begin(), edge_candidates.end()); SearchWorkspace workspace(edge_candidates); const int mincut_begin = static_cast(edge_candidates.size()); for (const auto& cut : generate_mincuts(input_seed)) workspace.add_cut(cut); const int candidate_count = static_cast(workspace.distances.size()); const int candidate_max = candidate_count - 1; auto random_candidate = [&]() { // 最小カット候補は経路候補より少ないので、選択確率を別に確保する。 if (candidate_count > mincut_begin && random_real() < MINCUT_PROB) { return randint(mincut_begin, candidate_max); } return randint(0, candidate_max); }; cerr << "mincut_candidates=" << candidate_count - mincut_begin << '\n'; vector edges, prev_edges, best_trees; edges.reserve(128); prev_edges.reserve(128); int est_sizes = 0; double current_score = numeric_limits::infinity(); long long best_score = numeric_limits::max(); auto contains_edge = [&](int edge) { return find(edges.begin(), edges.end(), edge) != edges.end(); }; while (true) { double elapsed = chrono::duration(chrono::steady_clock::now() - start_time).count(); if (elapsed > ROUTE_TIME_LIMIT) break; if (candidate_max < 0) break; double temperature = START_TEMP + (END_TEMP - START_TEMP) * elapsed / ROUTE_TIME_LIMIT; int x = random_candidate(); double r = random_real(); if (contains_edge(x)) continue; prev_edges = edges; int prev_est_sizes = est_sizes; double move_r = est_sizes >= 90 ? random_real() : 0.0; if (est_sizes >= 90 && move_r >= 1.0 - SHUFFLE_PROB) { shuffle(edges.begin(), edges.end(), rng); } else if (est_sizes >= 90 && move_r >= 1.0 - SHUFFLE_PROB - TWO_ADD_TWO_REMOVE_PROB) { if (edges.size() < 2 || candidate_count - static_cast(edges.size()) < 2) continue; int z = random_candidate(); while (z == x || contains_edge(z)) z = random_candidate(); int y1 = randint(0, static_cast(edges.size()) - 1); int y2 = randint(0, static_cast(edges.size()) - 2); if (y2 >= y1) ++y2; if (y1 < y2) swap(y1, y2); for (int y : {y1, y2}) { est_sizes -= workspace.distances[edges[y]]; edges.erase(edges.begin() + y); } for (int new_edge : {x, z}) { edges.push_back(new_edge); est_sizes += workspace.distances[new_edge]; } } else if (est_sizes <= 90) { est_sizes += workspace.distances[x]; edges.push_back(x); } else if (move_r < TWO_ADD_ONE_REMOVE_PROB && est_sizes < 105) { if (candidate_count - static_cast(edges.size()) < 2) continue; int z = random_candidate(); while (z == x || contains_edge(z)) z = random_candidate(); int y = randint(0, static_cast(edges.size()) - 1); est_sizes -= workspace.distances[edges[y]]; edges.erase(edges.begin() + y); for (int new_edge : {x, z}) { edges.push_back(new_edge); est_sizes += workspace.distances[new_edge]; } } else if (TWO_ADD_ONE_REMOVE_PROB <= move_r && move_r < TWO_ADD_ONE_REMOVE_PROB + TWO_REMOVE_ONE_ADD_PROB) { if (edges.size() < 2) continue; int y1 = randint(0, static_cast(edges.size()) - 1); int y2 = randint(0, static_cast(edges.size()) - 2); if (y2 >= y1) ++y2; if (y1 < y2) swap(y1, y2); for (int y : {y1, y2}) { est_sizes -= workspace.distances[edges[y]]; edges.erase(edges.begin() + y); } edges.push_back(x); est_sizes += workspace.distances[x]; } else { int y = randint(0, static_cast(edges.size()) - 1); est_sizes -= workspace.distances[edges[y]]; edges[y] = x; est_sizes += workspace.distances[x]; } const auto& trees = workspace.build(edges); long long new_score = workspace.evaluate(trees); if (new_score <= current_score || r < exp((current_score - new_score) / temperature)) { current_score = static_cast(new_score); } else { edges.swap(prev_edges); est_sizes = prev_est_sizes; } if (new_score < best_score) { best_trees = trees; best_score = new_score; } } refine_trees(best_trees, start_time); cout << best_trees.size() << '\n'; for (int v : best_trees) { auto [x, y] = to_2d(v); cout << x - 1 << ' ' << y - 1 << '\n'; } return 0; }