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

• 問題: No.5025 Forest Fire
• ユーザー: Shun_PI
• 提出日時: 2026-10-24 17:28:54
• 言語: C++23(gcc16) (gcc 16.2.0 + boost 1.92.0 + ACL)
• 結果: AC
• 実行時間: 1,737 ms / 2,000 ms
• コード長: 12,602 bytes
• コンパイル時間: 5,576 ms
• コンパイル使用メモリ: 397,144 KB
• 実行使用メモリ: 9,924 KB
• スコア: 181,411,617
• 最終ジャッジ日時: 2026-10-24 17:39:00
• 合計ジャッジ時間: 146,052 ms
• ジャッジサーバーID （参考情報）: judge6_0 / <nil>

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

ソースコード

#include <string>
#pragma GCC target("avx2")
#pragma GCC optimize("O3")
#pragma GCC optimize("unroll-loops")
#include <bits/stdc++.h>
using namespace std;
// MT19937 のパラメータ・出力列を保ち、状態の格納幅だけを 32 bit にする。
using SamplingRng = mersenne_twister_engine<uint32_t, 32, 624, 397, 31,
    0x9908b0dfU, 11, 0xffffffffU, 7, 0x9d2c5680U, 15, 0xefc60000U, 18, 1812433253U>;
using lint = long long int;
using P = pair<int, int>;
using PL = pair<lint, lint>;
#define FOR(i, begin, end) for(int i=(begin),i##_end_=(end);i<i##_end_;i++)
#define IFOR(i, begin, end) for(int i=(end)-1,i##_begin_=(begin);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<lint> RH_B = {1532834020, 1388622299};
vector<lint> RH_M = {2147482409, 2147478017};
constexpr int INF = 2147483647;
void yes(bool expr) {cout << (expr ? "Yes" : "No") << "\n";}
template<class T>void chmax(T &a, const T &b) { if (a<b) a=b; }
template<class T>void chmin(T &a, const T &b) { if (b<a) a=b; }

struct GridGraph {
    int n;
    int cells;
    vector<array<int, 4>> 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<int16_t, CAPACITY> positions;
    array<int16_t, CAPACITY> 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<int16_t>(padded);
                flat_index[padded] = static_cast<int16_t>(flat);
            }
        }
    }

    void accumulate(vector<long long> &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<int, 4> offsets;
    array<int16_t, CAPACITY> flat_index, low, subtree_size, separated_size, component;
    array<int, CAPACITY> separated_squares;
    array<unsigned char, CAPACITY> articulation;
    int timer = 0, component_count = 0;

    void dfs(int v, int parent) {
        const int own_order = ++timer;
        order[v] = static_cast<int16_t>(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<int16_t>(lowest);
        subtree_size[v] = static_cast<int16_t>(size);
        separated_size[v] = static_cast<int16_t>(separated);
        separated_squares[v] = squares;
        articulation[v] = parent == -1 ? parts >= 2 : parts >= 1;
        component[component_count++] = static_cast<int16_t>(v);
    }
};


// 提出時のサンプリング回数。実行時は --samples X でも変更できる。
constexpr int DEFAULT_SAMPLE_COUNT = 1000;
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<unsigned char> &dry) {
    vector<unsigned char> visited(graph.cells, false);
    vector<int> 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<long long>(component.size());
        total += size * size;
    }
    return total;
}

vector<int> select_best_cuts(const vector<long long> &decrease,
                            const vector<unsigned char> &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<int>(available.size()); ++v) {
            if (available[v] && (best == -1 || decrease[v] > decrease[best])) best = v;
        }
        if (best == -1) return {};
        return {best};
    }
    vector<int> candidates;
    candidates.reserve(available.size());
    for (int v = 0; v < static_cast<int>(available.size()); ++v) {
        if (available[v]) candidates.push_back(v);
    }
    const int count = min(max_count, static_cast<int>(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<unsigned char> blocked(available.size(), false);
    vector<int> selected;
    selected.reserve(count);
    while (!candidates.empty() && static_cast<int>(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<string> board(n);
    for (auto &row : board) cin >> row;

    GridGraph graph(n);
    vector<unsigned char> dryness_by_cell(graph.cells);
    vector<unsigned char> available(graph.cells);
    for (int v = 0; v < graph.cells; ++v) {
        dryness_by_cell[v] = static_cast<unsigned char>(board[v / n][v % n] - '0');
        available[v] = dryness_by_cell[v] != 0;
    }

    // 乱数生成器は一度だけ初期化し、各 step で新しい乱数を取得する。
    SamplingRng rng(0);
    uniform_int_distribution<int> random_digit(0, 9);
    uniform_int_distribution<int> 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<int> answer;
    vector<long long> 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<int, 10> dry_cutoff;
        for (int d = 0; d <= 9; ++d) {
            // R * 分母 < d * 分子 と、R < ceil(d * 分子 / 分母) は同値。
            dry_cutoff[d] = static_cast<int>((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<int>(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<unsigned char> 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;
}