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

• 問題: No.5025 Forest Fire
• ユーザー: nmnmnmnm
• 提出日時: 2026-10-24 17:28:03
• 言語: C++23(gcc16) (gcc 16.2.0 + boost 1.92.0 + ACL)
• 結果: AC
• 実行時間: 1,911 ms / 2,000 ms
• コード長: 20,891 bytes
• コンパイル時間: 4,410 ms
• コンパイル使用メモリ: 373,472 KB
• 実行使用メモリ: 9,924 KB
• スコア: 152,025,515
• 最終ジャッジ日時: 2026-10-24 17:39:11
• 合計ジャッジ時間: 161,109 ms
• ジャッジサーバーID （参考情報）: judge3_0 / <nil>

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

ソースコード

// https://github.com/nmnmnmnm04/ahc-template

#include <bits/stdc++.h>
#include <cassert>

using namespace std;

using ll = long long;
using ld = long double;
using u64 = unsigned long long;
using u128 = unsigned __int128;

#define rep(i, n) for (ll i = 0; i < (n); ++i)
#define drep(i, n) for (ll i = (n) - 1; i >= 0; --i)
#define nfor(i, s, n) for (ll i = (s); i < (n); ++i)
#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define endl "\n"
#define nall(v) (v).begin(), (v).end()
#define rall(v) (v).rbegin(), (v).rend()
#define pcnt(x) __builtin_popcountll(x)
#define next_p(v) next_permutation((v).begin(), (v).end())

template <class T, class U>
bool chmin(T& current, const U& candidate) {
    if (current <= candidate) return false;
    current = candidate;
    return true;
}

template <class T, class U>
bool chmax(T& current, const U& candidate) {
    if (current >= candidate) return false;
    current = candidate;
    return true;
}

const ll INF = 1LL << 61;

// 先頭4方向は U, D, L, R、後ろ4方向は斜め
const int di[8] = {-1, 1, 0, 0, -1, -1, 1, 1};
const int dj[8] = {0, 0, -1, 1, -1, 1, -1, 1};
const char dc[4] = {'U', 'D', 'L', 'R'};

bool in_grid(ll i, ll j, ll h, ll w) {
    return 0 <= i && i < h && 0 <= j && j < w;
}

#ifndef AHC_TEMPLATE_DEBUG
#define AHC_TEMPLATE_DEBUG 0
#endif
constexpr bool DEBUG = AHC_TEMPLATE_DEBUG;

class Timer {
public:
    using Clock = chrono::steady_clock;

    explicit Timer(double time_limit_seconds)
        : started_at_(Clock::now()), time_limit_seconds_(time_limit_seconds) {}

    void reset() { started_at_ = Clock::now(); }

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

    bool expired() const { return elapsed() >= time_limit_seconds_; }

    double time_limit() const { return time_limit_seconds_; }

    double remaining() const {
        return max(0.0, time_limit_seconds_ - elapsed());
    }

    double progress() const {
        if (time_limit_seconds_ <= 0.0) return 1.0;
        return min(1.0, max(0.0, elapsed() / time_limit_seconds_));
    }

private:
    Clock::time_point started_at_;
    double time_limit_seconds_;
};

// SplitMix64、高速かつシードを固定するとローカル実験を再現できる
class Random {
public:
    explicit Random(u64 seed = 0x123456789abcdef0ULL) : state_(seed) {}

    u64 next_u64() {
        u64 value = (state_ += 0x9e3779b97f4a7c15ULL);
        value = (value ^ (value >> 30)) * 0xbf58476d1ce4e5b9ULL;
        value = (value ^ (value >> 27)) * 0x94d049bb133111ebULL;
        return value ^ (value >> 31);
    }

    // [0, upper)
    int uniform_int(int upper) {
        assert(upper > 0);
        const u128 product = u128(next_u64()) * upper;
        return int(product >> 64);
    }

    // [lower, upper)
    int uniform_int(int lower, int upper) {
        assert(lower < upper);
        return lower + uniform_int(upper - lower);
    }

    // [0, 1)
    double uniform01() {
        return double(next_u64() >> 11) / 9007199254740992.0;
    }

private:
    u64 state_;
};

constexpr double TIME_LIMIT = 1.92;
Timer timer(TIME_LIMIT);
Random rng(0x123456789abcdef0ULL);

struct Forest {
    int n;
    int k;
    vector<string> d;
    vector<string> original_d;
    vector<int> danger;
    int total_start_weight = 0;
    vector<int> cut_cells;
    vector<vector<int>> cut_lines;

    Forest(vector<string> dryness, int budget)
        : n(int(dryness.size())), k(budget), d(move(dryness)), original_d(d),
          danger(n * n, 0) {
        rep(i, n) rep(j, n) {
            const int id = int(i * n + j);
            const int value = d[i][j] - '0';
            danger[id] = value * value;
            rep(direction, 4) {
                const int ni = int(i) + di[direction];
                const int nj = int(j) + dj[direction];
                if (!in_grid(ni, nj, n, n)) continue;
                const int neighbor = d[ni][nj] - '0';
                danger[id] += neighbor * neighbor;
            }
            if (value > 0) total_start_weight += danger[id];
        }
    }

    bool has_tree(int id) const {
        return 0 <= id && id < n * n && d[id / n][id % n] != '0';
    }

    // 抽選対象は未伐採の木のみ。重みは自身と4近傍の乾燥度の二乗和。
    int choose_start(Random& random) const {
        if (int(cut_cells.size()) >= k || total_start_weight == 0) return -1;
        int ticket = random.uniform_int(total_start_weight);
        int id = 0;
        for (const string& row : d) {
            for (const char value : row) {
                if (value != '0') {
                    if (ticket < danger[id]) return id;
                    ticket -= danger[id];
                }
                ++id;
            }
        }
        assert(false);
        return -1;
    }

    bool cut_tree(int id) {
        if (int(cut_cells.size()) >= k || !has_tree(id)) return false;
        const int i = id / n;
        const int j = id % n;
        const int value = d[i][j] - '0';
        const int removed = value * value;

        total_start_weight -= danger[id];
        d[i][j] = '0';
        danger[id] -= removed;
        rep(direction, 4) {
            const int ni = i + di[direction];
            const int nj = j + dj[direction];
            if (!in_grid(ni, nj, n, n)) continue;
            const int neighbor_id = ni * n + nj;
            danger[neighbor_id] -= removed;
            if (has_tree(neighbor_id)) total_start_weight -= removed;
        }
        cut_cells.push_back(id);
        return true;
    }

    void apply_cut_line(vector<int> line) {
        assert(!line.empty());
        for (const int id : line) {
            const bool cut = cut_tree(id);
            assert(cut);
            (void)cut;
        }
        cut_lines.push_back(move(line));
    }

    // 指定列の木を元に戻し、危険度・抽選重みを差分更新する。
    void remove_cut_line(int index) {
        assert(0 <= index && index < int(cut_lines.size()));
        const vector<int> removed_line = move(cut_lines[index]);
        cut_lines.erase(cut_lines.begin() + index);
        for (const int id : removed_line) {
            const int i = id / n;
            const int j = id % n;
            assert(d[i][j] == '0' && original_d[i][j] != '0');
            const int value = original_d[i][j] - '0';
            const int restored = value * value;
            d[i][j] = original_d[i][j];
            danger[id] += restored;
            total_start_weight += danger[id];
            rep(direction, 4) {
                const int ni = i + di[direction];
                const int nj = j + dj[direction];
                if (!in_grid(ni, nj, n, n)) continue;
                const int neighbor_id = ni * n + nj;
                danger[neighbor_id] += restored;
                if (has_tree(neighbor_id)) total_start_weight += restored;
            }
        }
        cut_cells.clear();
        for (const auto& line : cut_lines) {
            cut_cells.insert(cut_cells.end(), line.begin(), line.end());
        }
    }
};

constexpr array<int, 4> CUT_LENGTHS = {8, 16, 24, 32};
constexpr int DRY_SIMULATION_COUNT = 3;
constexpr double INITIAL_TEMPERATURE_RATIO = 0.01;
constexpr double FINAL_TEMPERATURE_RATIO = 0.00001;

// 方向は左上、右上、左下、右下の4種類。長さには開始地点も含める。
vector<int> make_cut_line(const Forest& forest, int start, int direction, int length) {
    assert(0 <= direction && direction < 4);
    vector<int> line;
    const int remaining = forest.k - int(forest.cut_cells.size());
    const int limit = min(length, remaining);
    if (limit <= 0) return line;
    line.reserve(limit);
    int current = start;
    while (int(line.size()) < limit && forest.has_tree(current)) {
        line.push_back(current);
        if (int(line.size()) == limit) break;
        const int ni = current / forest.n + di[direction + 4];
        const int nj = current % forest.n + dj[direction + 4];
        if (!in_grid(ni, nj, forest.n, forest.n)) break;
        current = ni * forest.n + nj;
    }
    return line;
}

vector<int> grow_cut_line(Forest& trial, int start, int direction, int length) {
    vector<int> line = make_cut_line(trial, start, direction, length);
    for (const int id : line) {
        const bool cut = trial.cut_tree(id);
        assert(cut);
        (void)cut;
    }
    return line;
}

// floor(random * 10 / 2^64) < d と同値な閾値。
// 木がないマスでも乱数を1回消費し、従来の乱数列を保つ。
constexpr array<u64, 10> DRY_THRESHOLDS = [] {
    array<u64, 10> thresholds{};
    for (int value = 1; value <= 9; ++value) {
        thresholds[value] = u64(((u128(value) << 64) + 9) / 10);
    }
    return thresholds;
}();

vector<unsigned char> sample_dryness(const Forest& forest, Random& random) {
    vector<unsigned char> dry(forest.n * forest.n);
    int id = 0;
    for (const string& row : forest.d) {
        for (const char value : row) {
            dry[id++] = random.next_u64() < DRY_THRESHOLDS[value - '0'];
        }
    }
    return dry;
}

// 外周を0で囲み、4近傍の境界判定と除算を省く。バッファは評価間で再利用。
struct SimulationWorkspace {
    int n = 0;
    int stride = 0;
    vector<unsigned char> dry;
    vector<int> queue;

    void resize(int side) {
        if (side == n) return;
        n = side;
        stride = n + 2;
        dry.assign(stride * stride, 0);
        queue.resize(n * n);
    }

    void sample(const Forest& forest, Random& random) {
        int offset = stride + 1;
        for (const string& row : forest.d) {
            for (const char value : row) {
                dry[offset++] = random.next_u64() < DRY_THRESHOLDS[value - '0'];
            }
            offset += 2;
        }
    }

    ll square_sum() {
        ll score = 0;
        unsigned char* active = dry.data();
        int* pending = queue.data();
        const int end = (n + 1) * stride;
        for (int start = stride + 1; start < end; ++start) {
            if (!active[start]) continue;
            int head = 0;
            int tail = 1;
            pending[0] = start;
            active[start] = 0;
            while (head < tail) {
                const int id = pending[head++];
                // 固定の4方向は展開して、反復ごとのループ制御を省く。
                const int up = id - stride;
                if (active[up]) {
                    active[up] = 0;
                    pending[tail++] = up;
                }
                const int down = id + stride;
                if (active[down]) {
                    active[down] = 0;
                    pending[tail++] = down;
                }
                const int left = id - 1;
                if (active[left]) {
                    active[left] = 0;
                    pending[tail++] = left;
                }
                const int right = id + 1;
                if (active[right]) {
                    active[right] = 0;
                    pending[tail++] = right;
                }
            }
            score += ll(tail) * tail;
        }
        return score;
    }
};

SimulationWorkspace& simulation_workspace(int n) {
    static thread_local SimulationWorkspace workspace;
    workspace.resize(n);
    return workspace;
}

// 乾燥した木の4近傍連結成分について、サイズの二乗和を返す。
ll component_square_sum(int n, const vector<unsigned char>& dry) {
    auto& workspace = simulation_workspace(n);
    for (int i = 0; i < n; ++i) {
        copy_n(dry.data() + i * n, n, workspace.dry.data() + (i + 1) * workspace.stride + 1);
    }
    return workspace.square_sum();
}

ll evaluate_forest(const Forest& forest, Random& random, const Timer* deadline = nullptr,
                   const vector<int>* additional_cuts = nullptr) {
    auto& workspace = simulation_workspace(forest.n);
    ll score = 0;
    rep(sample, DRY_SIMULATION_COUNT) {
        if (deadline != nullptr && deadline->expired()) return INF;
        workspace.sample(forest, random);
        if (additional_cuts != nullptr) {
            for (const int id : *additional_cuts) {
                const int padded_id = (id / forest.n + 1) * workspace.stride + id % forest.n + 1;
                workspace.dry[padded_id] = 0;
            }
        }
        score += workspace.square_sum();
    }
    // 全候補で共通の係数1/(DRY_SIMULATION_COUNT*N*N)は比較時に省略。
    return score;
}

struct ScoredLine {
    vector<int> cells;
    ll score = INF;
};

ScoredLine choose_cut_line(const Forest& forest, int start, Random& random,
                          const Timer* deadline = nullptr, int fixed_length = 0) {
    ScoredLine best;
    rep(direction, 4) {
        const int length_count = fixed_length > 0 ? 1 : int(CUT_LENGTHS.size());
        rep(length_index, length_count) {
            const int length = fixed_length > 0 ? fixed_length : CUT_LENGTHS[length_index];
            if (deadline != nullptr && deadline->expired()) return best;
            vector<int> line = make_cut_line(forest, start, int(direction), length);
            const ll score = evaluate_forest(forest, random, deadline, &line);
            // 同点では先に調べた候補を保持する。
            if (score < best.score) {
                best.score = score;
                best.cells = move(line);
            }
        }
    }
    return best;
}

void solve(Forest& forest, Random& random) {
    while (true) {
        const int start = forest.choose_start(random);
        if (start == -1) break;
        ScoredLine best = choose_cut_line(forest, start, random);
        forest.apply_cut_line(move(best.cells));
    }
}

array<int, 2> sample_split_lengths(int length, Random& random) {
    assert(length >= 10);
    const int center = length / 2;
    const int radius = length / 10 + 1;
    return {random.uniform_int(center - radius, center + radius + 1),
            random.uniform_int(center - radius, center + radius + 1)};
}

// 開始地点を危険度の重み付き抽選で選び、最良の列を設置する。
// fixed_length == 0 のときは初期解と同じ4種類の長さを試す。
ll place_best_line(Forest& forest, Random& random, const Timer* deadline = nullptr,
                   int fixed_length = 0) {
    if (deadline != nullptr && deadline->expired()) return INF;
    const int start = forest.choose_start(random);
    if (start == -1) return INF;
    ScoredLine best = choose_cut_line(forest, start, random, deadline, fixed_length);
    if (deadline != nullptr && deadline->expired()) return INF;
    if (best.cells.empty()) return INF;
    forest.apply_cut_line(move(best.cells));
    return best.score;
}

struct ExponentialTemperature {
    double initial_temperature;
    double log_cooling_ratio;

    explicit ExponentialTemperature(ll initial_score)
        : initial_temperature(double(initial_score) * INITIAL_TEMPERATURE_RATIO),
          log_cooling_ratio(log(FINAL_TEMPERATURE_RATIO / INITIAL_TEMPERATURE_RATIO)) {}

    double at(double progress) const {
        if (initial_temperature <= 0.0) return 0.0;
        progress = clamp(progress, 0.0, 1.0);
        return initial_temperature * exp(log_cooling_ratio * progress);
    }
};

bool accept_annealing_move(ll current_score, ll candidate_score, double temperature,
                           Random& random) {
    if (candidate_score == INF) return false;
    if (candidate_score <= current_score) return true;
    if (temperature <= 0.0) return false;
    const double delta = double(candidate_score - current_score);
    return random.uniform01() < exp(-delta / temperature);
}

// 既存の伐採位置から伸ばせる方向が残っているかを調べる。
bool has_cut_extension(const Forest& forest) {
    for (const int anchor : forest.cut_cells) {
        const int i = anchor / forest.n;
        const int j = anchor % forest.n;
        for (int direction = 4; direction < 8; ++direction) {
            const int ni = i + di[direction];
            const int nj = j + dj[direction];
            if (in_grid(ni, nj, forest.n, forest.n) && forest.has_tree(ni * forest.n + nj)) {
                return true;
            }
        }
    }
    return false;
}

int fill_remaining_cuts(Forest& forest, Random& random, const Timer* deadline = nullptr) {
    const int initial_count = int(forest.cut_cells.size());
    while (int(forest.cut_cells.size()) < forest.k && !forest.cut_cells.empty()) {
        if (deadline != nullptr && deadline->expired()) break;
        const int anchor = forest.cut_cells[random.uniform_int(int(forest.cut_cells.size()))];
        const int direction = random.uniform_int(4);
        const int ni = anchor / forest.n + di[direction + 4];
        const int nj = anchor % forest.n + dj[direction + 4];
        if (!in_grid(ni, nj, forest.n, forest.n) || !forest.has_tree(ni * forest.n + nj)) {
            if (!has_cut_extension(forest)) break;
            continue;
        }
        const int remaining = forest.k - int(forest.cut_cells.size());
        forest.apply_cut_line(make_cut_line(forest, ni * forest.n + nj, direction, remaining));
    }
    return int(forest.cut_cells.size()) - initial_count;
}

bool try_relocate_line(Forest& forest, int index, ll& current_score, Random& random,
                       const Timer* deadline = nullptr, double temperature = 0.0) {
    if (deadline != nullptr && deadline->expired()) return false;
    const int removed_length = int(forest.cut_lines[index].size());
    Forest base = forest;
    base.remove_cut_line(index);

    Forest candidate = base;
    ll candidate_score = place_best_line(candidate, random, deadline);
    if (removed_length >= 10) {
        if (deadline != nullptr && deadline->expired()) return false;
        const auto lengths = sample_split_lengths(removed_length, random);
        Forest split_candidate = base;
        const ll first_score = place_best_line(split_candidate, random, deadline, lengths[0]);
        if (first_score == INF) return false;
        const ll split_score = place_best_line(split_candidate, random, deadline, lengths[1]);
        if (split_score < candidate_score) {
            candidate_score = split_score;
            candidate = move(split_candidate);
        }
    }
    if (deadline != nullptr && deadline->expired()) return false;
    if (fill_remaining_cuts(candidate, random, deadline) > 0) {
        candidate_score = evaluate_forest(candidate, random, deadline);
    }
    if (deadline != nullptr && deadline->expired()) return false;
    if (!accept_annealing_move(current_score, candidate_score, temperature, random)) return false;
    forest = move(candidate);
    current_score = candidate_score;
    return true;
}

ll simulated_annealing(Forest& forest, Random& random, const Timer& deadline) {
    if (forest.cut_lines.empty() || deadline.expired()) return INF;
    ll current_score = evaluate_forest(forest, random, &deadline);
    if (current_score == INF) return INF;
    const ll initial_score = current_score;
    ll best_score = current_score;
    Forest best_forest = forest;
    const double started_at = deadline.elapsed();
    const double duration = deadline.time_limit() - started_at;
    if (duration <= 0.0) return best_score;
    const double inverse_duration = 1.0 / duration;
    const ExponentialTemperature temperature_schedule(initial_score);
    u64 iterations = 0;
    u64 accepted = 0;
    while (true) {
        const double now = deadline.elapsed();
        if (now >= deadline.time_limit()) break;
        ++iterations;
        const double progress = (now - started_at) * inverse_duration;
        const double temperature = temperature_schedule.at(progress);
        const int index = random.uniform_int(int(forest.cut_lines.size()));
        if (try_relocate_line(forest, index, current_score, random, &deadline, temperature)) {
            ++accepted;
            if (current_score < best_score) {
                best_score = current_score;
                best_forest = forest;
            }
        }
    }
    forest = move(best_forest);
    if constexpr (DEBUG) {
        cerr << "[sa] iterations=" << iterations << " accepted=" << accepted
             << " elapsed=" << deadline.elapsed() << '\n';
    }
    return best_score;
}

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

    int n, k;
    if (!(cin >> n >> k)) return 0;
    vector<string> d(n);
    for (string& row : d) cin >> row;
    Forest forest(move(d), k);

    solve(forest, rng);
    simulated_annealing(forest, rng, timer);
    cout << forest.cut_cells.size() << '\n';
    for (const int id : forest.cut_cells) {
        cout << id / n << ' ' << id % n << '\n';
    }
    return 0;
}
#endif