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

• 問題: No.5025 Forest Fire
• ユーザー: mech_39
• 提出日時: 2026-10-19 23:24:17
• 言語: C++23(gcc16) (gcc 16.2.0 + boost 1.92.0 + ACL)
• 結果: AC
• 実行時間: 1,841 ms / 2,000 ms
• コード長: 14,248 bytes
• コンパイル時間: 4,052 ms
• コンパイル使用メモリ: 383,260 KB
• 実行使用メモリ: 9,924 KB
• スコア: 183,879,289
• 最終ジャッジ日時: 2026-10-24 14:35:24
• 合計ジャッジ時間: 153,896 ms
• ジャッジサーバーID （参考情報）: judge4_0 / <nil>

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

ソースコード

#include <bits/stdc++.h>
#include <cassert>
using namespace std;

struct ForestRng {
    uint64_t x;

    explicit ForestRng(uint64_t s = 0x6a09e667f3bcc909ULL) : x(s) {}

    uint64_t operator()() {
        x ^= x << 7;
        x ^= x >> 9;
        return x;
    }

    int integer(int n) {
        return int((*this)() % uint64_t(n));
    }

    double real() {
        return double((*this)() >> 11) * 0x1.0p-53;
    }
};

struct ForestSamples {
    int N, V, maxSamples;
    vector<int> d, vertices;
    vector<array<int, 4>> adj;
    vector<uint8_t> degree, dry;

    vector<int> order, low, sub, sepSum, sepSq, component;
    const uint8_t *dryNow = nullptr;
    const uint8_t *cutNow = nullptr;
    int timer = 0;

    ForestSamples(
        const vector<string>& grid,
        int sampleCount = 512,
        uint64_t seed = 6523487656ULL
    )
        : N((int)grid.size()),
          V(N * N),
          maxSamples(sampleCount),
          d(V),
          adj(V),
          degree(V),
          dry(size_t(V) * sampleCount),
          order(V),
          low(V),
          sub(V),
          sepSum(V),
          sepSq(V),
          component(V) {

        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                int u = i * N + j;
                d[u] = grid[i][j] - '0';
                if (d[u]) vertices.push_back(u);
            }
        }

        for (int u : vertices) {
            int i = u / N, j = u % N;
            if (i && d[u - N])
                adj[u][degree[u]++] = u - N;
            if (j && d[u - 1])
                adj[u][degree[u]++] = u - 1;
            if (j + 1 < N && d[u + 1])
                adj[u][degree[u]++] = u + 1;
            if (i + 1 < N && d[u + N])
                adj[u][degree[u]++] = u + N;
        }

        // 各20サンプルで，各頂点を正確に2*d回乾燥させる．
        ForestRng rng(seed);
        const int block = 20;
        array<int, 20> idx;

        for (int base = 0; base < sampleCount; base += block) {
            int len = min(block, sampleCount - base);

            for (int u : vertices) {
                iota(idx.begin(), idx.end(), 0);

                int quota = d[u] * len / 10;
                if (d[u] * len % 10 &&
                    rng.integer(10) < d[u] * len % 10) {
                    ++quota;
                }

                for (int t = 0; t < quota; t++) {
                    swap(idx[t], idx[t + rng.integer(len - t)]);
                    dry[size_t(base + idx[t]) * V + u] = 1;
                }
            }
        }
    }

    void dfs(int u, int parent) {
        order[u] = low[u] = ++timer;
        component[timer - 1] = u;

        int sz = 1, ss = 0, sq = 0;

        for (int k = 0; k < degree[u]; k++) {
            int v = adj[u][k];

            if (v == parent || cutNow[v] || !dryNow[v])
                continue;

            if (!order[v]) {
                dfs(v, u);
                sz += sub[v];
                low[u] = min(low[u], low[v]);

                if (low[v] >= order[u]) {
                    ss += sub[v];
                    sq += sub[v] * sub[v];
                }
            } else {
                low[u] = min(low[u], order[v]);
            }
        }

        sub[u] = sz;
        sepSum[u] = ss;
        sepSq[u] = sq;
    }

    // gainsを指定すると，全頂点の削除利得も計算する．
    int64_t evaluate(
        const vector<uint8_t>& cuts,
        int samples,
        vector<int64_t>* gains = nullptr,
        int offset = 0
    ) {
        assert(samples + offset <= maxSamples);

        cutNow = cuts.data();
        if (gains) fill(gains->begin(), gains->end(), 0);

        int64_t total = 0;

        for (int s = offset; s < offset + samples; s++) {
            dryNow = dry.data() + size_t(s) * V;
            fill(order.begin(), order.end(), 0);
            timer = 0;

            if (gains) {
                for (int u : vertices) {
                    if (!dryNow[u] || cutNow[u] || order[u])
                        continue;

                    int begin = timer;
                    dfs(u, -1);

                    int sz = timer - begin;
                    total += sz * sz;

                    for (int z = begin; z < timer; z++) {
                        int v = component[z];
                        int rem = sz - 1 - sepSum[v];

                        (*gains)[v] +=
                            int64_t(sz) * sz
                            - sepSq[v]
                            - rem * rem;
                    }
                }
            } else {
                // 利得不要なら，通常の連結成分探索だけを行う．
                for (int u : vertices) {
                    if (!dryNow[u] || cutNow[u] || order[u])
                        continue;

                    int head = timer, begin = timer;
                    component[timer++] = u;
                    order[u] = 1;

                    while (head < timer) {
                        int v = component[head++];

                        for (int k = 0; k < degree[v]; k++) {
                            int w = adj[v][k];

                            if (!order[w] && dryNow[w] && !cutNow[w]) {
                                order[w] = 1;
                                component[timer++] = w;
                            }
                        }
                    }

                    int sz = timer - begin;
                    total += sz * sz;
                }
            }
        }

        return total;
    }
};

// 全て伐採した状態から，木を復元する構築法．
struct ReverseBuilder {
    ForestSamples& m;
    ForestRng rng;
    int S, V, offset;

    vector<int> parent, size;
    vector<uint8_t> cut;

    ReverseBuilder(
        ForestSamples& model,
        int samples,
        uint64_t seed,
        int off = 0
    )
        : m(model),
          rng(seed),
          S(samples),
          V(model.V),
          offset(off),
          parent(samples * V, -1),
          size(samples * V),
          cut(V) {

        for (int u : m.vertices) cut[u] = 1;
    }

    int root(int s, int u) {
        int off = s * V;

        while (parent[off + u] != u) {
            parent[off + u] = parent[off + parent[off + u]];
            u = parent[off + u];
        }

        return u;
    }

    // u自身は常に乾燥すると仮定して周囲との合併コストを求め，
    // 最後に乾燥度を掛ける．共通係数1/(10*S)は省略する．
    double cost(int u) {
        int64_t total = 0;

        for (int s = 0; s < S; s++) {
            const uint8_t* dry =
                m.dry.data() + size_t(s + offset) * V;

            int roots[4], count = 0;
            int merged = 1, oldSquares = 0;

            for (int k = 0; k < m.degree[u]; k++) {
                int v = m.adj[u][k];
                if (cut[v] || !dry[v]) continue;

                int r = root(s, v);
                bool duplicate = false;

                for (int j = 0; j < count; j++) {
                    if (roots[j] == r) duplicate = true;
                }

                if (duplicate) continue;

                roots[count++] = r;
                int z = size[s * V + r];
                merged += z;
                oldSquares += z * z;
            }

            total += merged * merged - oldSquares;
        }

        return double(total) * m.d[u];
    }

    void restore(int u) {
        cut[u] = 0;

        for (int s = 0; s < S; s++) {
            const uint8_t* dry =
                m.dry.data() + size_t(s + offset) * V;

            if (!dry[u]) continue;

            int off = s * V;
            parent[off + u] = u;
            size[off + u] = 1;

            for (int k = 0; k < m.degree[u]; k++) {
                int v = m.adj[u][k];
                if (cut[v] || !dry[v]) continue;

                int a = root(s, u), b = root(s, v);
                if (a == b) continue;

                if (size[off + a] < size[off + b])
                    swap(a, b);

                parent[off + b] = a;
                size[off + a] += size[off + b];
            }
        }
    }

    vector<uint8_t> run(int K) {
        vector<int> remaining = m.vertices;

        while ((int)remaining.size() > K) {
            int chosen = -1;
            double best = 1e100;

            bool exhaustive = remaining.size() <= 300;
            int count = exhaustive ? (int)remaining.size() : 40;

            for (int k = 0; k < count; k++) {
                int j = exhaustive
                    ? k
                    : rng.integer((int)remaining.size());

                double value =
                    cost(remaining[j]) * (0.999 + 0.002 * rng.real());

                if (value < best) {
                    best = value;
                    chosen = j;
                }
            }

            restore(remaining[chosen]);
            remaining[chosen] = remaining.back();
            remaining.pop_back();
        }

        return cut;
    }
};

struct Solver {
    using Clock = chrono::steady_clock;

    Clock::time_point start;
    double limit;
    ForestSamples model;
    ForestRng rng;
    int K;

    vector<vector<uint8_t>> archive;

    Solver(const vector<string>& grid, int k, double seconds)
        : start(Clock::now()),
          limit(seconds),
          model(grid, 1800),
          rng(987654321),
          K(min(k, (int)model.vertices.size())) {}

    double elapsed() const {
        return chrono::duration<double>(Clock::now() - start).count();
    }

    int64_t validation(const vector<uint8_t>& cut, int samples = 200) {
        return model.evaluate(cut, samples, nullptr, 1000);
    }

    vector<uint8_t> refine(vector<uint8_t> cut, double deadline) {
        vector<int> removed;
        for (int u : model.vertices) {
            if (cut[u]) removed.push_back(u);
        }

        vector<int64_t> gain(model.V);
        int S = 20;

        int64_t score = model.evaluate(cut, S);
        int64_t bestScore = score;
        auto best = cut;

        double begin = elapsed();
        double duration = max(1e-6, deadline - begin);

        while (elapsed() < deadline && !removed.empty()) {
            double progress = (elapsed() - begin) / duration;

            int target =
                progress < .22 ? 20 :
                progress < .52 ? 40 :
                progress < .78 ? 80 : 160;

            if (target != S) {
                archive.push_back(best);
                cut = best;
                removed.clear();

                for (int u : model.vertices) {
                    if (cut[u]) removed.push_back(u);
                }

                S = target;
                score = bestScore = model.evaluate(cut, S);
            }

            int index = rng.integer((int)removed.size());
            int add = removed[index];
            cut[add] = 0;

            int64_t restored = model.evaluate(cut, S, &gain);

            int del = -1;
            for (int u : model.vertices) {
                if (cut[u] || u == add) continue;

                if (del < 0 ||
                    gain[u] > gain[del] ||
                    (gain[u] == gain[del] && (rng() & 7) == 0)) {
                    del = u;
                }
            }

            if (del < 0) {
                cut[add] = 1;
                break;
            }

            int64_t next = restored - gain[del];

            double phase = fmod(progress * 4, 1.0);
            double temperature = score * .0015 * pow(.02, phase);
            if (progress > .88) temperature *= .05;

            if (next <= score ||
                rng.real() <
                    exp(double(score - next) / max(1e-9, temperature))) {

                cut[del] = 1;
                removed[index] = del;
                score = next;

                if (score < bestScore) {
                    bestScore = score;
                    best = cut;
                }
            } else {
                cut[add] = 1;
            }
        }

        archive.push_back(best);
        return best;
    }

    vector<uint8_t> solve() {
        if (K == (int)model.vertices.size()) {
            vector<uint8_t> cut(model.V);
            for (int u : model.vertices) cut[u] = 1;
            return cut;
        }

        if (K == 0) return vector<uint8_t>(model.V);

        vector<pair<int64_t, vector<uint8_t>>> candidates;
        int attempts = 0;

        do {
            ReverseBuilder builder(
                model,
                20,
                987654321ULL + 55434ULL * attempts,
                0
            );

            auto cut = builder.run(K);
            candidates.emplace_back(validation(cut), move(cut));
            ++attempts;

        } while (elapsed() < limit * .25);

        sort(candidates.begin(), candidates.end());

        auto best = candidates[0].second;
        archive.push_back(best);

        best = refine(best, limit * .91);

        vector<pair<int64_t, int>> ranked;
        for (int i = 0; i < (int)archive.size(); i++) {
            ranked.emplace_back(validation(archive[i]), i);
        }

        sort(ranked.begin(), ranked.end());

        int64_t bestScore = LLONG_MAX;
        for (int j = 0; j < min(3, (int)ranked.size()); j++) {
            int idx = ranked[j].second;
            int64_t value = validation(archive[idx], 800);

            if (value < bestScore) {
                bestScore = value;
                best = archive[idx];
            }
        }

        cerr << "starts=" << attempts
             << " elapsed=" << elapsed()
             << " seconds\n";

        return best;
    }
};

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    int N, K;
    if (!(cin >> N >> K)) return 0;

    vector<string> grid(N);
    for (auto& row : grid) cin >> row;

    double seconds = 1.85;
    if (const char* s = getenv("TIME_LIMIT")) {
        seconds = max(.1, atof(s));
    }

    Solver solver(grid, K, seconds);
    auto cut = solver.solve();

    cout << accumulate(cut.begin(), cut.end(), 0) << '\n';

    for (int u = 0; u < N * N; u++) {
        if (cut[u]) {
            cout << u / N << ' ' << u % N << '\n';
        }
    }
}