#1195113 (C++23) No.5025 Forest Fire

• 問題: No.5025 Forest Fire
• ユーザー: Piiiii
• 提出日時: 2026-10-24 17:08:04
• 言語: C++23 (gcc 15.3.0 + boost 1.92.0 + ACL)
• 結果: AC
• 実行時間: 1,967 ms / 2,000 ms
• コード長: 18,263 bytes
• コンパイル時間: 4,361 ms
• コンパイル使用メモリ: 377,496 KB
• 実行使用メモリ: 9,928 KB
• スコア: 146,600,133
• 最終ジャッジ日時: 2026-10-24 17:10:49
• 合計ジャッジ時間: 160,035 ms
• ジャッジサーバーID （参考情報）: judge6_0 / <nil>

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

ソースコード

#pragma GCC optimize("O3")
#include <bits/stdc++.h>
using namespace std;

// tools.zip/src/lib.rs の採点処理を、先頭10サンプルについて移植する。
// 採用する操作は、指定された「伐採位置を1つ、別の乾燥木へ移す」のみ。
// スコアが厳密に改善した操作だけを採用する山登り法。

constexpr int SAMPLES = 10;
constexpr double SEARCH_TIME_LIMIT_MS = 1900.0;
constexpr int STALL_TRIALS = 10000;
constexpr int FINAL_SAMPLES = 1000;
constexpr uint64_t GAMMA = 0x9E3779B97F4A7C15ULL;

uint64_t splitmix64(uint64_t &state) {
    state += GAMMA;
    uint64_t z = state;
    z = (z ^ (z >> 30)) * 0xBF58476D1CE4E5B9ULL;
    z = (z ^ (z >> 27)) * 0x94D049BB133111EBULL;
    return z ^ (z >> 31);
}

struct Board {
    int n, k, v;
    vector<uint8_t> d;
    vector<uint8_t> cut;
    vector<array<int, 4>> neighbors;
    vector<uint8_t> degree;

    Board(int n_, int k_, vector<uint8_t> d_)
        : n(n_), k(k_), v(n * n), d(move(d_)), cut(v, 0),
          neighbors(v), degree(v, 0) {
        for (int r = 0; r < n; ++r) {
            for (int c = 0; c < n; ++c) {
                const int x = r * n + c;
                if (r > 0) neighbors[x][degree[x]++] = x - n;
                if (r + 1 < n) neighbors[x][degree[x]++] = x + n;
                if (c > 0) neighbors[x][degree[x]++] = x - 1;
                if (c + 1 < n) neighbors[x][degree[x]++] = x + 1;
            }
        }
    }

    // Rust の sample_seed(input) と同一。
    uint64_t sample_seed() const {
        uint64_t h = 1469598103934665603ULL;
        for (uint8_t x : d) {
            h ^= uint64_t('0' + x);
            h *= 1099511628211ULL;
        }
        return h ^ (uint64_t(n) * 1000003ULL + uint64_t(k) * 998244353ULL);
    }

    // Rust の sample_dry(input, sample_index) と同一。
    vector<uint8_t> generate_dry(uint64_t sample_index) const {
        uint64_t state = sample_seed()
                       + uint64_t(sample_index) * uint64_t((v + 1) / 2) * GAMMA;
        vector<uint8_t> dry(v, 0);
        for (int x = 0; x < v; x += 2) {
            const uint64_t r = splitmix64(state);
            const uint64_t lo_threshold = uint64_t(double(d[x]) * (4294967296.0 / 10.0));
            dry[x] = uint8_t((r & 0xFFFFFFFFULL) < lo_threshold);
            if (x + 1 < v) {
                const uint64_t hi_threshold = uint64_t(double(d[x + 1]) * (4294967296.0 / 10.0));
                dry[x + 1] = uint8_t((r >> 32) < hi_threshold);
            }
        }
        return dry;
    }
};

// 伐採なしの1サンプルにおける、連結成分のサイズの二乗和。
long long score_without_cuts(const Board &b, const vector<uint8_t> &dry) {
    vector<uint8_t> seen(b.v, 0);
    vector<int> queue;
    queue.reserve(b.v);
    long long total = 0;
    for (int x = 0; x < b.v; ++x) {
        if (!dry[x] || seen[x]) continue;
        queue.clear();
        queue.push_back(x);
        seen[x] = 1;
        for (size_t h = 0; h < queue.size(); ++h) {
            int u = queue[h];
            for (int j = 0; j < b.degree[u]; ++j) {
                int w = b.neighbors[u][j];
                if (dry[w] && !seen[w]) {
                    seen[w] = 1;
                    queue.push_back(w);
                }
            }
        }
        total += 1LL * queue.size() * queue.size();
    }
    return total;
}

// 現在の連結成分ごとに、各マスを取り除いたときの分割サイズを保持する。
// 候補の評価では連結成分の探索を行わず、採用された場合だけ再構築する。
// 評価関数と、交換候補の選択・採用条件は変更しない。
struct SampleComponents {
    const Board &b;
    vector<uint8_t> dry;
    vector<int> component, component_size;
    vector<int> tin, tout, low, parent, subtree_size;
    vector<int> rest_size;
    vector<array<int, 4>> separating_child;
    vector<uint8_t> separating_count;
    vector<long long> deletion_delta;
    long long q = 0;
    int timer = 0;

    SampleComponents(const Board &board, vector<uint8_t> dry_)
        : b(board), dry(move(dry_)),
          component(b.v, -1), tin(b.v, -1), tout(b.v), low(b.v),
          parent(b.v, -1), subtree_size(b.v), rest_size(b.v),
          separating_child(b.v), separating_count(b.v, 0), deletion_delta(b.v) {
        rebuild();
    }

    void rebuild() {
        fill(tin.begin(), tin.end(), -1);
        fill(component.begin(), component.end(), -1);
        fill(parent.begin(), parent.end(), -1);
        component_size.clear();
        timer = 0;
        q = 0;
        // DFSの木と到達可能な辺を使って、各頂点を取り除いた場合の分割を調べる。
        auto dfs = [&](auto &&self, int u, int id) -> void {
            component[u] = id;
            tin[u] = low[u] = timer++;
            subtree_size[u] = 1;
            for (int j = 0; j < b.degree[u]; ++j) {
                int w = b.neighbors[u][j];
                if (!dry[w] || b.cut[w] || w == parent[u]) continue;
                if (tin[w] < 0) {
                    parent[w] = u;
                    self(self, w, id);
                    subtree_size[u] += subtree_size[w];
                    low[u] = min(low[u], low[w]);
                } else {
                    low[u] = min(low[u], tin[w]);
                }
            }
            tout[u] = timer;
        };
        for (int u = 0; u < b.v; ++u) {
            if (!dry[u] || b.cut[u] || tin[u] >= 0) continue;
            int id = (int)component_size.size();
            dfs(dfs, u, id);
            int sz = subtree_size[u];
            component_size.push_back(sz);
            q += 1LL * sz * sz;
        }
        for (int u = 0; u < b.v; ++u) {
            if (component[u] < 0) continue;
            const int sz = component_size[component[u]];
            int apart = 0;
            long long after_sum = 0;
            int count = 0;
            for (int j = 0; j < b.degree[u]; ++j) {
                int w = b.neighbors[u][j];
                if (parent[w] == u && low[w] >= tin[u]) {
                    separating_child[u][count++] = w;
                    apart += subtree_size[w];
                    after_sum += 1LL * subtree_size[w] * subtree_size[w];
                }
            }
            separating_count[u] = (uint8_t)count;
            rest_size[u] = sz - 1 - apart;
            after_sum += 1LL * rest_size[u] * rest_size[u];
            deletion_delta[u] = after_sum - 1LL * sz * sz;
        }
    }

    // メモリ上の成分情報だけで、交換による正確な評価値の差分を返す。
    long long evaluate_swap(int add, int remove) const {
        const int old_id = component[remove]; // sample上で非乾燥なら -1
        long long delta = old_id < 0 ? 0 : deletion_delta[remove];
        if (!dry[add]) return delta;

        int used_keys[4];
        int used_count = 0;
        long long joined_size = 1, before_join = 0;
        for (int j = 0; j < b.degree[add]; ++j) {
            int w = b.neighbors[add][j];
            if (w == remove || !dry[w] || b.cut[w]) continue;
            int id = component[w];
            int key = id, sz = component_size[id];
            if (id == old_id) {
                // removeによって切り離されるDFSの子の部分木に含まれるかを調べる。
                key = old_id;
                sz = rest_size[remove];
                for (int t = 0; t < separating_count[remove]; ++t) {
                    int child = separating_child[remove][t];
                    if (tin[child] <= tin[w] && tin[w] < tout[child]) {
                        key = -1 - t;
                        sz = subtree_size[child];
                        break;
                    }
                }
            }
            bool unique = true;
            for (int t = 0; t < used_count; ++t) {
                if (used_keys[t] == key) { unique = false; break; }
            }
            if (!unique) continue;
            used_keys[used_count++] = key;
            joined_size += sz;
            before_join += 1LL * sz * sz;
        }
        return delta + joined_size * joined_size - before_join;
    }
};

long long official_score_10(long long q0, long long q, int v) {
    // tools.zip の Simulation::e, e0, score と同じ計算。
    const double e0 = (double)q0 / SAMPLES / v;
    const double e  = (double)q / SAMPLES / v;
    return llround(1e6 * (1.0 + log2(e0 / max(e, 1e-12))));
}


// 候補の比較に使う1000サンプルは、tools.zip と同じ生成方法で得る。
// 各行の乾燥マスを64bit整数に記録して、全候補について使い回す。
// 連結成分のサイズ二乗和を、各行内の連続する木をまとめて計算する。
struct FinalEvaluator {
    const Board &b;
    vector<uint64_t> sample_rows;
    vector<int> parent, size;

    explicit FinalEvaluator(const Board &board)
        : b(board), sample_rows((size_t)FINAL_SAMPLES * b.n),
          parent(b.v), size(b.v) {
        const uint64_t seed = b.sample_seed();
        const uint64_t advance = uint64_t((b.v + 1) / 2) * GAMMA;
        vector<uint64_t> thresholds(b.v);
        for (int x = 0; x < b.v; ++x) {
            thresholds[x] = uint64_t(double(b.d[x]) * (4294967296.0 / 10.0));
        }
        if (b.n % 2 == 0) {
            // 偶数幅では64bit乱数の上位・下位32bitが同じ行に収まる。
            for (int s = 0; s < FINAL_SAMPLES; ++s) {
                uint64_t state = seed + uint64_t(s) * advance;
                uint64_t *rows = sample_rows.data() + (size_t)s * b.n;
                for (int i = 0; i < b.n; ++i) {
                    uint64_t mask = 0;
                    const uint64_t *thr = thresholds.data() + i * b.n;
                    for (int j = 0; j < b.n; j += 2) {
                        const uint64_t rnd = splitmix64(state);
                        mask |= uint64_t((rnd & 0xffffffffULL) < thr[j]) << j;
                        mask |= uint64_t((rnd >> 32) < thr[j + 1]) << (j + 1);
                    }
                    rows[i] = mask;
                }
            }
        } else {
            // 奇数幅では乱数のペアが行境界をまたぐため、元と同じ順序で処理。
            for (int s = 0; s < FINAL_SAMPLES; ++s) {
                uint64_t state = seed + uint64_t(s) * advance;
                uint64_t *rows = sample_rows.data() + (size_t)s * b.n;
                for (int x = 0; x < b.v; x += 2) {
                    const uint64_t rnd = splitmix64(state);
                    if ((rnd & 0xffffffffULL) < thresholds[x]) {
                        rows[x / b.n] |= 1ULL << (x % b.n);
                    }
                    if (x + 1 < b.v && (rnd >> 32) < thresholds[x + 1]) {
                        rows[(x + 1) / b.n] |= 1ULL << ((x + 1) % b.n);
                    }
                }
            }
        }
    }

    inline int find(int x) {
        while (parent[x] != x) {
            parent[x] = parent[parent[x]];
            x = parent[x];
        }
        return x;
    }

    inline long long unite(int a, int c) {
        a = find(a);
        c = find(c);
        if (a == c) return 0;
        if (size[a] < size[c]) swap(a, c);
        long long delta = 2LL * size[a] * size[c];
        parent[c] = a;
        size[a] += size[c];
        return delta;
    }

    long long evaluate(const vector<int> &cuts) {
        uint64_t cut_rows[64] = {};
        for (int v : cuts) cut_rows[v / b.n] |= 1ULL << (v % b.n);

        long long total_q = 0;
        for (int s = 0; s < FINAL_SAMPLES; ++s) {
            const uint64_t *row = sample_rows.data() + (size_t)s * b.n;
            uint64_t prev_mask = 0, prev_starts = 0;
            long long q = 0;
            for (int i = 0; i < b.n; ++i) {
                const uint64_t mask = row[i] & ~cut_rows[i];
                const uint64_t starts = mask & ~(mask << 1);
                uint64_t x = starts;
                while (x) {
                    const int j = __builtin_ctzll(x);
                    x &= x - 1;
                    const uint64_t suffix = mask >> j;
                    const int len = (~suffix) ? __builtin_ctzll(~suffix) : 64 - j;
                    const int id = i * b.n + j;
                    parent[id] = id;
                    size[id] = len;
                    q += 1LL * len * len;
                }
                uint64_t both_starts = (prev_mask & mask) & ~((prev_mask & mask) << 1);
                while (both_starts) {
                    const int j = __builtin_ctzll(both_starts);
                    both_starts &= both_starts - 1;
                    const uint64_t prefix = ~0ULL >> (63 - j);
                    const int prev_j = 63 - __builtin_clzll(prev_starts & prefix);
                    const int now_j = 63 - __builtin_clzll(starts & prefix);
                    q += unite((i - 1) * b.n + prev_j, i * b.n + now_j);
                }
                prev_mask = mask;
                prev_starts = starts;
            }
            total_q += q;
        }
        return total_q;
    }
};

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

    int n, k;
    if (!(cin >> n >> k)) return 0;
    vector<uint8_t> d(n * n);
    for (int i = 0; i < n; ++i) {
        string row;
        cin >> row;
        for (int j = 0; j < n; ++j) d[i * n + j] = uint8_t(row[j] - '0');
    }
    Board board(n, k, move(d));

    vector<int> candidate_cells;
    candidate_cells.reserve(board.v);
    for (int x = 0; x < board.v; ++x) {
        if (board.d[x] >= 6) candidate_cells.push_back(x);
    }
    mt19937_64 rng((uint64_t)chrono::steady_clock::now().time_since_epoch().count());
    shuffle(candidate_cells.begin(), candidate_cells.end(), rng);
    int m = min(k, (int)candidate_cells.size());
    vector<int> cut_cells(candidate_cells.begin(), candidate_cells.begin() + m);
    vector<int> uncut_cells(candidate_cells.begin() + m, candidate_cells.end());
    for (int x : cut_cells) board.cut[x] = 1;

    // 現在の伐採位置を引き継ぎ、各10サンプルの山登り終了時の最良解を保存。
    vector<vector<int>> saved_cuts;
    const auto start = chrono::steady_clock::now();
#ifdef HC_STATS
    uint64_t trials = 0, accepted = 0, stalls = 0;
#endif

    if (!cut_cells.empty() && !uncut_cells.empty()) {
        vector<SampleComponents> sample_states;
        sample_states.reserve(SAMPLES);
        long long q0_total = 0, q_total = 0;
        uint64_t next_sample_index = 0;

        auto load_samples = [&]() {
            sample_states.clear();
            q0_total = q_total = 0;
            for (int s = 0; s < SAMPLES; ++s) {
                vector<uint8_t> dry = board.generate_dry(next_sample_index++);
                q0_total += score_without_cuts(board, dry);
                sample_states.emplace_back(board, move(dry));
                q_total += sample_states.back().q;
            }
        };
        load_samples();
        long long current_score = official_score_10(q0_total, q_total, board.v);
        int unsuccessful = 0;

        uniform_int_distribution<int> cut_dist(0, (int)cut_cells.size() - 1);
        uniform_int_distribution<int> uncut_dist(0, (int)uncut_cells.size() - 1);
        while (chrono::duration<double, milli>(
                   chrono::steady_clock::now() - start).count() < SEARCH_TIME_LIMIT_MS) {
            const int ci = cut_dist(rng);
            const int ui = uncut_dist(rng);
            const int add = cut_cells[ci];
            const int remove = uncut_cells[ui];
            long long delta_total = 0;
            for (int s = 0; s < SAMPLES; ++s) {
                delta_total += sample_states[s].evaluate_swap(add, remove);
            }
            const long long candidate_score = official_score_10(
                q0_total, q_total + delta_total, board.v
            );
#ifdef HC_STATS
            ++trials;
#endif
            if (candidate_score > current_score) {
                board.cut[add] = 0;
                board.cut[remove] = 1;
                for (int s = 0; s < SAMPLES; ++s) {
                    if (sample_states[s].dry[add] || sample_states[s].dry[remove]) {
                        sample_states[s].rebuild();
                    }
                }
                swap(cut_cells[ci], uncut_cells[ui]);
                q_total += delta_total;
                current_score = candidate_score;
                unsuccessful = 0;
#ifdef HC_STATS
                ++accepted;
#endif
            } else if (++unsuccessful == STALL_TRIALS) {
                saved_cuts.push_back(cut_cells);
                load_samples(); // 次の10サンプルを取り直し、現在の解から続ける。
                current_score = official_score_10(q0_total, q_total, board.v);
                unsuccessful = 0;
#ifdef HC_STATS
                ++stalls;
#endif
            }
        }
    }
    saved_cuts.push_back(cut_cells); // 最後の10サンプルでの最良解も保存する。

    // 1000サンプルをすべての保存解に共通で用い、期待燃焼数が最小の解を選ぶ。
    // 同じ解が複数回保存された場合、評価結果を再利用する。
    // 交換で配列の順序だけが変わる場合もあるので、各集合をソートして同一判定する。
    for (auto &cuts : saved_cuts) sort(cuts.begin(), cuts.end());
    sort(saved_cuts.begin(), saved_cuts.end());
    saved_cuts.erase(unique(saved_cuts.begin(), saved_cuts.end()), saved_cuts.end());

    FinalEvaluator evaluator(board);
    long long best_q = LLONG_MAX;
    vector<int> best_cuts;
    for (const auto &cuts : saved_cuts) {
        long long q = evaluator.evaluate(cuts);
        if (q < best_q) {
            best_q = q;
            best_cuts = cuts;
        }
    }

#ifdef HC_STATS
    cerr << "trials=" << trials << " accepted=" << accepted
         << " resamples=" << stalls << " saved=" << saved_cuts.size()
         << " final_q=" << best_q << " elapsed_ms="
         << chrono::duration<double, milli>(chrono::steady_clock::now() - start).count()
         << '\n';
#endif
    cout << best_cuts.size() << '\n';
    for (int x : best_cuts) cout << x / n << ' ' << x % n << '\n';
    return 0;
}