#1195271 (C++23(gcc16)) No.5025 Forest Fire

• 問題: No.5025 Forest Fire
• ユーザー: prussian_coder
• 提出日時: 2026-10-24 22:11:48
• 言語: C++23(gcc16) (gcc 16.2.0 + boost 1.92.0 + ACL)
• 結果: AC
• 実行時間: 1,793 ms / 2,000 ms
• コード長: 19,036 bytes
• コンパイル時間: 4,372 ms
• コンパイル使用メモリ: 298,968 KB
• 実行使用メモリ: 11,000 KB
• スコア: 187,327,357
• 最終ジャッジ日時: 2026-10-24 22:14:21
• 合計ジャッジ時間: 152,954 ms
• ジャッジサーバーID （参考情報）: judge7_0 / <nil>

ファイルパターン	結果
other	AC * 50

ソースコード

#include <algorithm>
#include <array>
#include <chrono>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <limits>
#include <random>
#include <string>
#include <utility>
#include <unordered_set>
#include <vector>

using namespace std;
using Connection = pair<int, int>;

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 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<vector<int>> board;
vector<int> danger_zones;

int to_1d(int i, int j) { return i * N2 + j; }
pair<int, int> 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<int> 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<int> 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 SearchWorkspace {
    vector<int> routes, route_offsets;
    vector<int> distances, danger_idx, trees, parents;
    vector<pair<int, int>> adjacent_dangers;
    vector<unsigned int> used;
    vector<uint64_t> removal_mask, cached_removal_mask;
    vector<int> removed_indices;
    long long cached_score = 0;
    bool score_cached = false;
    unsigned int tree_stamp = 0;

    explicit SearchWorkspace(const vector<Connection>& 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<int>(routes.size()));
            distances.push_back(connection_dist(v1, v2));
        }
        for (int i = 0; i < static_cast<int>(danger_zones.size()); ++i) {
            danger_idx[danger_zones[i]] = i;
        }
        for (int i = 0; i < static_cast<int>(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<int>& build(const vector<int>& 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<int>(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<int>& selected_trees) {
        constexpr int REMOVED = numeric_limits<int>::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<int> vertices;
    vector<array<int, 4>> neighbors;
    vector<vector<unsigned char>> dry;
    vector<unsigned char> cut;
    vector<int> tin, low, subtree, component, component_sizes, order, separated;
    vector<long long> separated_sq, gain, restore_cost;
    int timer = 0;

    SampleEvaluator() : cells(N2 * N2), neighbors(cells),
        dry(SAMPLES, vector<unsigned char>(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<unsigned char>& 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<int>(component_sizes.size());
                int begin = static_cast<int>(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<int>(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<int, 4> 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;
    }
};

void refine_trees(vector<int>& trees, chrono::steady_clock::time_point start) {
    SampleEvaluator evaluator;
    for (int v : trees) evaluator.cut[v] = 1;
    auto expired = [&]() {
        return chrono::duration<double>(chrono::steady_clock::now() - start).count() >= TIME_LIMIT;
    };
    long long score = evaluator.analyze();
    vector<int> best_trees = trees;
    long long best_score = score;
    vector<unsigned char> tried(evaluator.cells, 0);
    int iterations = 0, improvements = 0;
    vector<long long> saved_restore_cost;
    while (!expired()) {
        bool changed = true;
        // 本数が足りないときは、期待焼失本数を最も減らす木から追加する。
        if (static_cast<int>(trees.size()) < K) {
            int add = -1;
            for (int v : evaluator.vertices) if (!evaluator.cut[v] &&
                (add < 0 || evaluator.gain[v] > evaluator.gain[add])) add = v;
            if (add < 0) break;
            evaluator.cut[add] = 1;
            trees.push_back(add);
            score -= evaluator.gain[add];
        } else {
            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 cost = evaluator.restore_cost[remove];
            saved_restore_cost = evaluator.restore_cost;
            evaluator.cut[remove] = 0;
            evaluator.analyze();
            int add = remove;
            for (int v : evaluator.vertices) if (!evaluator.cut[v] &&
                evaluator.gain[v] > evaluator.gain[add]) add = v;
            // 復元前後で同じサンプルを使い、厳密に改善する交換だけ採用する。
            if (evaluator.gain[add] > cost) {
                evaluator.cut[add] = 1;
                *find(trees.begin(), trees.end(), remove) = add;
                score += cost - evaluator.gain[add];
                fill(tried.begin(), tried.end(), 0);
                ++improvements;
            } else {
                evaluator.cut[remove] = 1;
                tried[remove] = 1;
                evaluator.restore_cost.swap(saved_restore_cost);
                changed = false;
            }
        }
        if (score < best_score) {
            best_score = score;
            best_trees = trees;
        }
        ++iterations;
        if (expired()) break;
        if (changed) evaluator.analyze();
    }
    trees = best_trees;
    cerr << "refine_iterations=" << iterations << " improvements=" << improvements
         << " sample_score=" << best_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<int>(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<uint32_t>(d)) * 16777619U;
    }
    mt19937 rng(input_seed);
    auto randint = [&](int lo, int hi) {
        return uniform_int_distribution<int>(lo, hi)(rng);
    };
    auto random_real = [&]() {
        return uniform_real_distribution<double>(0.0, 1.0)(rng);
    };
    vector<int> 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<uint64_t> candidate_set;
    candidate_set.reserve(CAND_SIZES_LIM + 1);
    vector<Connection> 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<int>(safe_zones.size()) - 1)];
        int v2 = safe_zones[randint(0, static_cast<int>(safe_zones.size()) - 1)];
        int dist = connection_dist(v1, v2);
        if (2 <= dist && dist <= 50) {
            uint64_t key = (uint64_t(static_cast<uint32_t>(v1)) << 32) |
                           static_cast<uint32_t>(v2);
            if (candidate_set.insert(key).second) edge_candidates.emplace_back(v1, v2);
        }
    }
    // 旧版の set と同じ候補順に戻し、乱数に対する選択結果を維持する。
    sort(edge_candidates.begin(), edge_candidates.end());
    const int candidate_max = static_cast<int>(edge_candidates.size()) - 1;
    SearchWorkspace workspace(edge_candidates);
    vector<int> edges, prev_edges, best_trees;
    edges.reserve(128);
    prev_edges.reserve(128);
    int est_sizes = 0;
    double current_score = numeric_limits<double>::infinity();
    long long best_score = numeric_limits<long long>::max();
    auto contains_edge = [&](int edge) {
        return find(edges.begin(), edges.end(), edge) != edges.end();
    };

    while (true) {
        double elapsed = chrono::duration<double>(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 = randint(0, candidate_max);
        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 || edge_candidates.size() - edges.size() < 2) continue;
            int z = randint(0, candidate_max);
            while (z == x || contains_edge(z)) z = randint(0, candidate_max);
            int y1 = randint(0, static_cast<int>(edges.size()) - 1);
            int y2 = randint(0, static_cast<int>(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 (edge_candidates.size() - edges.size() < 2) continue;
            int z = randint(0, candidate_max);
            while (z == x || contains_edge(z)) z = randint(0, candidate_max);
            int y = randint(0, static_cast<int>(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<int>(edges.size()) - 1);
            int y2 = randint(0, static_cast<int>(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<int>(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<double>(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;
}