#1195093 (Rust) No.5025 Forest Fire

• 問題: No.5025 Forest Fire
• ユーザー: g4np0n_kyopro
• 提出日時: 2026-10-24 16:53:55
• 言語: Rust (1.99.0 + proconio + num + itertools + ACL)
• 結果: AC
• 実行時間: 1,920 ms / 2,000 ms
• コード長: 37,039 bytes
• コンパイル時間: 2,071 ms
• コンパイル使用メモリ: 207,012 KB
• 実行使用メモリ: 9,936 KB
• スコア: 132,967,772
• 最終ジャッジ日時: 2026-10-24 16:56:32
• 合計ジャッジ時間: 157,428 ms
• ジャッジサーバーID （参考情報）: judge2_0 / <nil>

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

ソースコード

use std::collections::hash_map::RandomState;
use std::hash::{BuildHasher, Hasher};
use std::io::{self, Read};
use std::sync::atomic::{AtomicU64, Ordering};
use std::sync::OnceLock;
use std::time::{Duration, Instant};

// 定数定義
pub const EXPECTATION_SAMPLES: usize = 50;
const ANNEALING_TIME_MS: u64 = 1950;
const START_TEMPERATURE: f64 = 10.0;
const END_TEMPERATURE: f64 = 0.1;
const START_PENALTY: f64 = 0.1;
const END_PENALTY: f64 = 1000.0;
// 領域の上書き、辺の上書き、ON/OFF辺の入れ替え。
const NEIGHBOR_WEIGHTS: [usize; 3] = [1, 4, 4];
const NEIGHBOR_NAMES: [&str; 3] = ["area", "edge", "edge_swap"];
const AREA_BOUNDARIES: [usize; 7] = [5, 11, 17, 23, 29, 36, 43];
const AREA_COUNT_PER_AXIS: usize = AREA_BOUNDARIES.len() + 1;
const NONE: usize = usize::MAX;

fn next_expectation_seed() -> u64 {
    // 起動ごとにランダムな初期値を用意し、呼び出しごとにシードを進める。
    static STATE: OnceLock<AtomicU64> = OnceLock::new();
    let state = STATE.get_or_init(|| AtomicU64::new(RandomState::new().build_hasher().finish()));
    let mut seed = state.fetch_add(0x9E3779B97F4A7C15, Ordering::Relaxed);
    splitmix64(&mut seed)
}

fn splitmix64(state: &mut u64) -> u64 {
    *state = state.wrapping_add(0x9E3779B97F4A7C15);
    let mut z = *state;
    z = (z ^ (z >> 30)).wrapping_mul(0xBF58476D1CE4E5B9);
    z = (z ^ (z >> 27)).wrapping_mul(0x94D049BB133111EB);
    z ^ (z >> 31)
}

/// 各乾燥連結成分の「全体の本数×対象の本数」の和を返す。
/// dryの外周はfalseとし、訪問したマスをfalseに変更する。
fn component_target_sum(
    dry: &mut [bool],
    targets: &[bool],
    width: usize,
    stack: &mut Vec<usize>,
) -> u64 {
    let mut sum = 0;
    for i in 1..width - 1 {
        for j in 1..width - 1 {
            let start = i * width + j;
            if !dry[start] {
                continue;
            }
            dry[start] = false;
            stack.push(start);
            let mut size = 0_u64;
            let mut target_size = 0_u64;
            while let Some(v) = stack.pop() {
                size += 1;
                target_size += u64::from(targets[v]);
                for next in [v - width, v + width, v - 1, v + 1] {
                    if dry[next] {
                        dry[next] = false;
                        stack.push(next);
                    }
                }
            }
            sum += size * target_size;
        }
    }
    sum
}

/// 期待焼失本数をモンテカルロ推定する。
/// drynessはmainと同じ、外周を0で埋めた(N+2)×(N+2)の盤面。
/// sample_countには通常EXPECTATION_SAMPLESを渡す。呼び出すたびにシードを変更する。
/// 各サンプルのΣ|C_k|²を、サンプル数とN²で割った値を返す。
/// 計算量はO(sample_count * N²)、追加メモリはO(N²)。
pub fn estimate_expected_burned(dryness: &[Vec<u32>], sample_count: usize) -> f64 {
    estimate_expected_burned_with_seed(dryness, None, sample_count, next_expectation_seed())
}

/// 指定集合に属する木の期待焼失本数を推定する。
/// treesは外周を含めない0始まりの座標(i, j)の集合。重複は1本として扱う。
/// 木のないマスや伐採済みのマスを含めても、その寄与は0になる。
/// 集合外の木も延焼を伝え、落雷位置は盤面全体から選ばれる。
/// 各サンプルでΣ(|C_k| * |C_k ∩ trees|) / N²を計算して平均する。
/// 盤面表現・サンプル数・シード更新はestimate_expected_burnedと同じ。
/// 現時点では全盤面をサンプリング・探索するため、計算量はO(sample_count * N²)。
pub fn estimate_expected_burned_for_trees(
    dryness: &[Vec<u32>],
    trees: &[(usize, usize)],
    sample_count: usize,
) -> f64 {
    estimate_expected_burned_with_seed(dryness, Some(trees), sample_count, next_expectation_seed())
}

// 固定シードの内部処理は、テストの再現性を保つために分離する。
fn estimate_expected_burned_with_seed(
    dryness: &[Vec<u32>],
    trees: Option<&[(usize, usize)]>,
    sample_count: usize,
    seed: u64,
) -> f64 {
    assert!(sample_count > 0, "サンプル数は1以上が必要です");
    let width = dryness.len();
    assert!(width >= 3, "盤面には1マス以上と外周が必要です");
    assert!(dryness.iter().all(|row| row.len() == width));
    let n = width - 2;
    let mut targets = vec![trees.is_none(); width * width];
    if let Some(trees) = trees {
        for &(i, j) in trees {
            assert!(i < n && j < n, "対象の座標は盤面内である必要があります");
            targets[(i + 1) * width + j + 1] = true;
        }
    }
    let mut thresholds = vec![0_u64; width * width];
    for i in 0..width {
        for j in 0..width {
            let d = dryness[i][j];
            assert!(d <= 9);
            if i == 0 || j == 0 || i == width - 1 || j == width - 1 {
                assert_eq!(d, 0, "外周は0で埋めてください");
            }
            // 配布ツールと同じ32ビット乱数による乾燥判定。
            thresholds[i * width + j] = (u64::from(d) << 32) / 10;
        }
    }

    let mut state = seed;
    let mut dry = vec![false; width * width];
    let mut stack = Vec::with_capacity(n * n);
    let mut total = 0.0;
    for _ in 0..sample_count {
        for i in 1..=n {
            for j in 1..=n {
                let v = i * width + j;
                let random = splitmix64(&mut state) as u32;
                dry[v] = u64::from(random) < thresholds[v];
            }
        }
        total += component_target_sum(&mut dry, &targets, width, &mut stack) as f64;
    }
    total / sample_count as f64 / (n * n) as f64
}

struct Board {
    n: usize,
    dryness: Vec<u32>,
    neighbors: Vec<[usize; 4]>,
}

impl Board {
    fn new(dryness: &[Vec<u32>]) -> Self {
        let n = dryness.len() - 2;
        let mut values = Vec::with_capacity(n * n);
        let mut neighbors = Vec::with_capacity(n * n);
        for i in 0..n {
            for j in 0..n {
                values.push(dryness[i + 1][j + 1]);
                neighbors.push([
                    if i > 0 { (i - 1) * n + j } else { NONE },
                    if i + 1 < n { (i + 1) * n + j } else { NONE },
                    if j > 0 { i * n + j - 1 } else { NONE },
                    if j + 1 < n { i * n + j + 1 } else { NONE },
                ]);
            }
        }
        Self {
            n,
            dryness: values,
            neighbors,
        }
    }
}

struct Layout {
    areas: Vec<Vec<usize>>,
    edges: Vec<Vec<usize>>,
    candidates: Vec<usize>,
}

impl Layout {
    fn new(board: &Board) -> Self {
        let count = AREA_COUNT_PER_AXIS;
        assert!(board.n > AREA_BOUNDARIES[count - 2] + 1);
        let starts: Vec<usize> = std::iter::once(0)
            .chain(AREA_BOUNDARIES.iter().map(|&line| line + 1))
            .collect();
        let ends: Vec<usize> = AREA_BOUNDARIES
            .iter()
            .copied()
            .chain(std::iter::once(board.n))
            .collect();

        let mut edges = Vec::new();
        // 横辺7×8本、縦辺8×7本。内部境界を盤面端まで延ばす。
        // 外周を囲む辺は作らず、内部境界同士の交点を除外する。
        // 共有辺は一度だけ登録。
        for &row in &AREA_BOUNDARIES {
            for col in 0..count {
                edges.push(
                    (starts[col]..ends[col])
                        .map(|j| row * board.n + j)
                        .filter(|&v| board.dryness[v] > 0)
                        .collect::<Vec<_>>(),
                );
            }
        }
        let vertical_offset = edges.len();
        for row in 0..count {
            for &col in &AREA_BOUNDARIES {
                edges.push(
                    (starts[row]..ends[row])
                        .map(|i| i * board.n + col)
                        .filter(|&v| board.dryness[v] > 0)
                        .collect::<Vec<_>>(),
                );
            }
        }
        let mut areas = Vec::new();
        for row in 0..count {
            for col in 0..count {
                let mut edge_ids = Vec::with_capacity(4);
                if row > 0 {
                    edge_ids.push((row - 1) * count + col);
                }
                if row + 1 < count {
                    edge_ids.push(row * count + col);
                }
                if col > 0 {
                    edge_ids.push(vertical_offset + row * (count - 1) + col - 1);
                }
                if col + 1 < count {
                    edge_ids.push(vertical_offset + row * (count - 1) + col);
                }
                areas.push(
                    edge_ids
                        .iter()
                        .flat_map(|&id| edges[id].iter().copied())
                        .collect(),
                );
            }
        }
        let candidates = edges.iter().flatten().copied().collect();
        Self {
            areas,
            edges,
            candidates,
        }
    }
}

// 乾燥状態によらず、現時点で木が存在するマスの4近傍連結成分を取り出す。
// 別成分間では延焼しないため、各成分の期待焼失本数を加算できる。
fn split_tree_sets(board: &Board, mut active: Vec<bool>) -> Vec<Vec<usize>> {
    let mut sets = Vec::new();
    for start in 0..active.len() {
        if !active[start] {
            continue;
        }
        active[start] = false;
        let mut cells = vec![start];
        let mut cursor = 0;
        while cursor < cells.len() {
            for &next in &board.neighbors[cells[cursor]] {
                if next != NONE && active[next] {
                    active[next] = false;
                    cells.push(next);
                }
            }
            cursor += 1;
        }
        sets.push(cells);
    }
    sets
}

// 他の木と接続していない集合の期待値。集合内だけサンプリング・探索する。
// 分母は集合サイズではなく、落雷先の総数である盤面全体のN²。
fn estimate_tree_set(
    board: &Board,
    cells: &[usize],
    sample_count: usize,
    seed: u64,
    deadline: Instant,
) -> Option<f64> {
    assert!(sample_count > 0);
    let mut index = vec![NONE; board.n * board.n];
    for (local, &v) in cells.iter().enumerate() {
        index[v] = local;
    }
    let neighbors: Vec<[usize; 4]> = cells
        .iter()
        .map(|&v| board.neighbors[v].map(|next| if next == NONE { NONE } else { index[next] }))
        .collect();
    let thresholds: Vec<u64> = cells
        .iter()
        .map(|&v| (u64::from(board.dryness[v]) << 32) / 10)
        .collect();
    let mut dry = vec![false; cells.len()];
    let mut stack = Vec::with_capacity(cells.len());
    let mut state = seed;
    let mut total = 0.0;
    for _ in 0..sample_count {
        if Instant::now() >= deadline {
            return None;
        }
        for (v, &threshold) in thresholds.iter().enumerate() {
            dry[v] = u64::from(splitmix64(&mut state) as u32) < threshold;
        }
        for start in 0..cells.len() {
            if !dry[start] {
                continue;
            }
            dry[start] = false;
            stack.push(start);
            let mut size = 0_u64;
            while let Some(v) = stack.pop() {
                size += 1;
                for &next in &neighbors[v] {
                    if next != NONE && dry[next] {
                        dry[next] = false;
                        stack.push(next);
                    }
                }
            }
            total += (size * size) as f64;
        }
    }
    Some(total / sample_count as f64 / (board.n * board.n) as f64)
}

struct TreeSet {
    cells: Vec<usize>,
    expected: f64,
}

fn evaluate_tree_sets(board: &Board, active: Vec<bool>, deadline: Instant) -> Option<Vec<TreeSet>> {
    split_tree_sets(board, active)
        .into_iter()
        .map(|cells| {
            let expected = estimate_tree_set(
                board,
                &cells,
                EXPECTATION_SAMPLES,
                next_expectation_seed(),
                deadline,
            )?;
            Some(TreeSet { cells, expected })
        })
        .collect()
}

struct ForestState {
    cuts: Vec<bool>,
    cut_count: usize,
    labels: Vec<usize>,
    sets: Vec<Option<TreeSet>>,
    free_ids: Vec<usize>,
    expected: f64,
}

struct Proposal {
    cuts: Vec<bool>,
    cut_count: usize,
    removed_ids: Vec<usize>,
    added_sets: Vec<TreeSet>,
    expected: f64,
}

enum ProposalResult {
    Unchanged,
    TimedOut,
    Changed(Proposal),
}

impl ForestState {
    fn new(board: &Board, candidates: &[usize], deadline: Instant) -> Option<Self> {
        let mut cuts = vec![false; board.n * board.n];
        for &v in candidates {
            cuts[v] = true;
        }
        let cut_count = cuts.iter().filter(|&&cut| cut).count();
        let active = board
            .dryness
            .iter()
            .enumerate()
            .map(|(v, &d)| d > 0 && !cuts[v])
            .collect();
        let sets = evaluate_tree_sets(board, active, deadline)?;
        let mut labels = vec![NONE; board.n * board.n];
        for (id, set) in sets.iter().enumerate() {
            for &v in &set.cells {
                labels[v] = id;
            }
        }
        let expected = sets.iter().map(|set| set.expected).sum();
        Some(Self {
            cuts,
            cut_count,
            labels,
            sets: sets.into_iter().map(Some).collect(),
            free_ids: Vec::new(),
            expected,
        })
    }

    fn propose(
        &self,
        board: &Board,
        selected: &[usize],
        on: bool,
        deadline: Instant,
    ) -> ProposalResult {
        self.propose_overwrites(board, &[(selected, on)], deadline)
    }

    // 複数辺の変更を一つの近傍として評価し、受理時にまとめて反映する。
    fn propose_overwrites(
        &self,
        board: &Board,
        updates: &[(&[usize], bool)],
        deadline: Instant,
    ) -> ProposalResult {
        let mut cuts = self.cuts.clone();
        for &(selected, on) in updates {
            for &v in selected {
                cuts[v] = on;
            }
        }
        let changed: Vec<usize> = cuts
            .iter()
            .enumerate()
            .filter_map(|(v, &cut)| (self.cuts[v] != cut).then_some(v))
            .collect();
        if changed.is_empty() {
            return ProposalResult::Unchanged;
        }
        let mut cut_count = self.cut_count;
        let mut affected = vec![false; self.sets.len()];
        for &v in &changed {
            if cuts[v] {
                cut_count += 1;
            } else {
                cut_count -= 1;
            }
            // 削除される木の旧集合と、復活する木に接続する旧集合を含める。
            for cell in std::iter::once(v).chain(board.neighbors[v].iter().copied()) {
                if cell != NONE && self.labels[cell] != NONE {
                    affected[self.labels[cell]] = true;
                }
            }
        }
        let removed_ids: Vec<usize> = (0..affected.len()).filter(|&id| affected[id]).collect();
        let mut active = vec![false; board.n * board.n];
        let mut old_expected = 0.0;
        for &id in &removed_ids {
            let set = self.sets[id].as_ref().unwrap();
            old_expected += set.expected;
            for &v in &set.cells {
                active[v] = !cuts[v];
            }
        }
        for &v in &changed {
            active[v] = !cuts[v];
        }
        let Some(added_sets) = evaluate_tree_sets(board, active, deadline) else {
            return ProposalResult::TimedOut;
        };
        let new_expected: f64 = added_sets.iter().map(|set| set.expected).sum();
        ProposalResult::Changed(Proposal {
            cuts,
            cut_count,
            removed_ids,
            added_sets,
            expected: self.expected - old_expected + new_expected,
        })
    }

    fn accept(&mut self, proposal: Proposal) {
        for id in proposal.removed_ids {
            let set = self.sets[id].take().unwrap();
            for v in set.cells {
                self.labels[v] = NONE;
            }
            self.free_ids.push(id);
        }
        for set in proposal.added_sets {
            let id = self.free_ids.pop().unwrap_or_else(|| {
                self.sets.push(None);
                self.sets.len() - 1
            });
            for &v in &set.cells {
                self.labels[v] = id;
            }
            self.sets[id] = Some(set);
        }
        self.cuts = proposal.cuts;
        self.cut_count = proposal.cut_count;
        self.expected = proposal.expected;
    }
}

#[derive(Default)]
struct NeighborStats {
    proposed: usize,
    evaluated: usize,
    accepted: usize,
    unchanged: usize,
    timed_out: usize,
}

fn annealing_parameters(progress: f64) -> (f64, f64) {
    let t = progress.clamp(0.0, 1.0);
    (
        START_TEMPERATURE + (END_TEMPERATURE - START_TEMPERATURE) * t,
        START_PENALTY * (END_PENALTY / START_PENALTY).powf(t),
    )
}

fn random_index(state: &mut u64, len: usize) -> usize {
    assert!(len > 0);
    // 剰余による偏りを避ける。
    let bound = len as u64;
    let threshold = bound.wrapping_neg() % bound;
    loop {
        let value = splitmix64(state);
        if value >= threshold {
            return (value % bound) as usize;
        }
    }
}

fn choose_edge_swap(layout: &Layout, cuts: &[bool], rng: &mut u64) -> Option<(usize, usize)> {
    let mut on_edges = Vec::new();
    let mut off_edges = Vec::new();
    for (id, edge) in layout.edges.iter().enumerate() {
        if edge.is_empty() {
            continue;
        }
        if edge.iter().all(|&v| cuts[v]) {
            on_edges.push(id);
        } else if edge.iter().all(|&v| !cuts[v]) {
            off_edges.push(id);
        }
    }
    if on_edges.is_empty() || off_edges.is_empty() {
        return None;
    }
    Some((
        on_edges[random_index(rng, on_edges.len())],
        off_edges[random_index(rng, off_edges.len())],
    ))
}

fn anneal(board: &Board, layout: &Layout, k: usize, started: Instant) -> Vec<usize> {
    let duration = Duration::from_millis(ANNEALING_TIME_MS);
    let deadline = started + duration;
    let mut state = ForestState::new(board, &layout.candidates, deadline)
        .expect("初期解の評価中に制限時間に到達しました");
    let mut rng = next_expectation_seed();
    let mut stats: [NeighborStats; 3] = std::array::from_fn(|_| NeighborStats::default());
    let total_weight: usize = NEIGHBOR_WEIGHTS.iter().sum();
    let (_, initial_penalty) =
        annealing_parameters(started.elapsed().as_secs_f64() / duration.as_secs_f64());
    // 各時点のペナルティ込み評価値の過去最小。探索状態の選択には使用しない。
    let mut best_score =
        state.expected + state.cut_count.saturating_sub(k) as f64 * initial_penalty;
    eprintln!("score\tbest_score\tneighbor\ton\tstatus\tcandidate_score\texpected_burned\tcuts\ttemperature\tpenalty");
    while Instant::now() < deadline && !layout.candidates.is_empty() {
        let draw = random_index(&mut rng, total_weight);
        let kind = if draw < NEIGHBOR_WEIGHTS[0] {
            0
        } else if draw < NEIGHBOR_WEIGHTS[0] + NEIGHBOR_WEIGHTS[1] {
            1
        } else {
            2
        };
        stats[kind].proposed += 1;
        let (result, on) = if kind == 2 {
            let result = if let Some((old_edge, new_edge)) =
                choose_edge_swap(layout, &state.cuts, &mut rng)
            {
                state.propose_overwrites(
                    board,
                    &[
                        (&layout.edges[old_edge], false),
                        (&layout.edges[new_edge], true),
                    ],
                    deadline,
                )
            } else {
                ProposalResult::Unchanged
            };
            (result, "NA")
        } else {
            let groups = if kind == 0 {
                &layout.areas
            } else {
                &layout.edges
            };
            let selected = &groups[random_index(&mut rng, groups.len())];
            let on = random_index(&mut rng, 2) == 1;
            (
                state.propose(board, selected, on, deadline),
                if on { "true" } else { "false" },
            )
        };
        let elapsed = started.elapsed();
        let (temperature, penalty) =
            annealing_parameters(elapsed.as_secs_f64() / duration.as_secs_f64());
        // 比較する双方の超過伐採ペナルティには同じ時刻の係数を使用。
        let old_score = state.expected + state.cut_count.saturating_sub(k) as f64 * penalty;
        let mut candidate_score = None;
        let mut timed_out = false;
        let status = match result {
            ProposalResult::Unchanged => {
                stats[kind].unchanged += 1;
                "unchanged"
            }
            ProposalResult::TimedOut => {
                stats[kind].timed_out += 1;
                timed_out = true;
                "timed_out"
            }
            ProposalResult::Changed(proposal) => {
                let new_score =
                    proposal.expected + proposal.cut_count.saturating_sub(k) as f64 * penalty;
                candidate_score = Some(new_score);
                if elapsed >= duration {
                    stats[kind].timed_out += 1;
                    timed_out = true;
                    "timed_out"
                } else {
                    stats[kind].evaluated += 1;
                    let delta = new_score - old_score;
                    let uniform = (splitmix64(&mut rng) >> 11) as f64 / ((1_u64 << 53) as f64);
                    let accepted = delta <= 0.0 || uniform < (-delta / temperature).exp();
                    if accepted {
                        stats[kind].accepted += 1;
                        state.accept(proposal);
                        "accepted"
                    } else {
                        "rejected"
                    }
                }
            }
        };
        let score = state.expected + state.cut_count.saturating_sub(k) as f64 * penalty;
        best_score = best_score.min(score);
        let name = NEIGHBOR_NAMES[kind];
        let candidate =
            candidate_score.map_or_else(|| "NA".to_owned(), |value| format!("{value:.6}"));
        eprintln!(
            "{score:.6}\t{best_score:.6}\t{name}\t{on}\t{status}\t{candidate}\t{:.6}\t{}\t{temperature:.6}\t{penalty:.6}",
            state.expected, state.cut_count
        );
        if timed_out {
            break;
        }
    }
    eprintln!("\nneighbor\tproposed\tevaluated\taccepted\tacceptance_rate\tunchanged\ttimed_out");
    for (name, stat) in NEIGHBOR_NAMES.iter().zip(&stats) {
        let rate = if stat.evaluated == 0 {
            0.0
        } else {
            stat.accepted as f64 / stat.evaluated as f64
        };
        eprintln!(
            "{name}\t{}\t{}\t{}\t{:.6}\t{}\t{}",
            stat.proposed, stat.evaluated, stat.accepted, rate, stat.unchanged, stat.timed_out
        );
    }
    let (_, penalty) =
        annealing_parameters(started.elapsed().as_secs_f64() / duration.as_secs_f64());
    eprintln!("\ncuts\tlimit\texpected_burned\tpenalty\tsamples");
    eprintln!(
        "{}\t{}\t{:.6}\t{:.6}\t{}",
        state.cut_count, k, state.expected, penalty, EXPECTATION_SAMPLES
    );
    assert!(
        state.cut_count <= k,
        "伐採本数が上限を超えています: {} > {}",
        state.cut_count,
        k
    );
    state
        .cuts
        .iter()
        .enumerate()
        .filter_map(|(v, &cut)| cut.then_some(v))
        .collect()
}

fn main() {
    // 入力の受け取りを始める前から計測する。
    let started = Instant::now();
    let mut input = String::new();
    io::stdin().read_to_string(&mut input).unwrap();
    let mut tokens = input.split_whitespace();
    let n: usize = tokens.next().unwrap().parse().unwrap();
    let k: usize = tokens.next().unwrap().parse().unwrap();
    // 外周を0で埋め、盤面外の乾燥度は合計に含めない。
    let mut dryness = vec![vec![0_u32; n + 2]; n + 2];
    for row in dryness.iter_mut().skip(1).take(n) {
        let digits = tokens.next().unwrap().as_bytes();
        for j in 0..n {
            row[j + 1] = (digits[j] - b'0') as u32;
        }
    }

    let board = Board::new(&dryness);
    let layout = Layout::new(&board);
    let cuts = anneal(&board, &layout, k, started);
    println!("{}", cuts.len());
    for v in cuts {
        println!("{} {}", v / n, v % n);
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    fn padded(rows: &[&[u32]]) -> Vec<Vec<u32>> {
        let n = rows.len();
        let mut board = vec![vec![0; n + 2]; n + 2];
        for (i, row) in rows.iter().enumerate() {
            board[i + 1][1..=n].copy_from_slice(row);
        }
        board
    }

    #[test]
    fn empty_board_has_zero_expectation() {
        assert_eq!(estimate_expected_burned(&padded(&[&[0]]), 10), 0.0);
    }

    #[test]
    fn components_use_four_neighbors() {
        let mut dry = vec![false; 25];
        // 隣接した2マスと、斜めに離れた1マス。二乗和は4+1。
        for (i, j) in [(1, 1), (1, 2), (2, 3)] {
            dry[i * 5 + j] = true;
        }
        assert_eq!(
            component_target_sum(&mut dry, &[true; 25], 5, &mut Vec::new()),
            5
        );
        assert!(dry.iter().all(|&cell| !cell));
    }

    #[test]
    fn sampling_matches_small_boards_with_known_expectations() {
        // 精度の検証に使う件数は、実行用の定数変更から独立させる。
        let samples = 40_000;
        let single = padded(&[&[7]]);
        let estimate = estimate_expected_burned_with_seed(&single, None, samples, 123);
        assert!((estimate - 0.7).abs() < 0.01);
        assert_eq!(
            estimate,
            estimate_expected_burned_with_seed(&single, None, samples, 123)
        );

        // 2×2盤面で隣接する乾燥度5の木が2本。
        // Σ|C_k|²の期待値は0.5*1 + 0.25*4 = 1.5、E=1.5/4。
        let adjacent = padded(&[&[5, 5], &[0, 0]]);
        let estimate = estimate_expected_burned_with_seed(&adjacent, None, samples, 123);
        assert!((estimate - 0.375).abs() < 0.01);
    }

    #[test]
    fn target_can_burn_from_lightning_outside_the_set() {
        // 乾燥度5の木が2本隣接。片方だけを対象とするとE=(0.25*1+0.25*2)/4。
        let board = padded(&[&[5, 5], &[0, 0]]);
        let estimate = estimate_expected_burned_with_seed(&board, Some(&[(0, 0)]), 40_000, 123);
        assert!((estimate - 0.1875).abs() < 0.01);
    }

    #[test]
    fn subset_expectations_partition_the_whole_on_shared_samples() {
        let board = padded(&[&[7, 8], &[3, 5]]);
        let evaluate = |trees| estimate_expected_burned_with_seed(&board, trees, 1000, 123);
        let whole = evaluate(None);
        assert_eq!(whole, evaluate(Some(&[(0, 0), (0, 1), (1, 0), (1, 1)])));
        let first = evaluate(Some(&[(0, 0), (1, 1)]));
        let second = evaluate(Some(&[(0, 1), (1, 0)]));
        assert!((whole - first - second).abs() < 1e-12);
        assert_eq!(evaluate(Some(&[(0, 0)])), evaluate(Some(&[(0, 0), (0, 0)])));
    }

    #[test]
    fn empty_and_removed_targets_have_zero_expectation() {
        let board = padded(&[&[9, 0], &[0, 9]]);
        assert_eq!(estimate_expected_burned_for_trees(&board, &[], 10), 0.0);
        assert_eq!(
            estimate_expected_burned_for_trees(&board, &[(0, 1)], 10),
            0.0
        );
    }

    #[test]
    fn successive_seeds_are_distinct() {
        let mut seeds = std::collections::HashSet::new();
        for _ in 0..10_000 {
            assert!(seeds.insert(next_expectation_seed()));
        }
    }

    #[test]
    #[should_panic(expected = "サンプル数は1以上")]
    fn zero_samples_are_rejected() {
        estimate_expected_burned(&padded(&[&[0]]), 0);
    }

    #[test]
    fn layout_has_shared_edges_and_excludes_corners_and_water() {
        let mut cells = vec![vec![1; 52]; 52];
        cells[6][2] = 0; // 盤面の(5, 1)は水。
        let board = Board::new(&cells);
        let layout = Layout::new(&board);
        assert_eq!(layout.areas.len(), 64);
        assert_eq!(layout.edges.len(), 112);
        let unique: std::collections::HashSet<_> = layout.candidates.iter().copied().collect();
        assert_eq!(unique.len(), layout.candidates.len());
        assert_eq!(unique.len(), 601); // 境界候補602マスから水1マスを除く。
        assert!(!unique.contains(&(5 * 50 + 1)));
        for i in 0..50 {
            for j in 0..50 {
                let row_boundary = [5, 11, 17, 23, 29, 36, 43].contains(&i);
                let col_boundary = [5, 11, 17, 23, 29, 36, 43].contains(&j);
                assert_eq!(
                    unique.contains(&(i * 50 + j)),
                    (row_boundary ^ col_boundary) && board.dryness[i * 50 + j] > 0
                );
            }
        }
        // 隣り合う領域の共有辺は、同じ5マスを操作する。
        let shared: Vec<_> = layout.areas[0]
            .iter()
            .filter(|v| layout.areas[1].contains(v))
            .copied()
            .collect();
        assert_eq!(shared.len(), 5);
        assert_eq!(layout.areas[63].len(), 12); // 右下6×6領域の内側2辺のみ。
        for group in layout.areas.iter().chain(&layout.edges) {
            assert!(group.iter().all(|&v| unique.contains(&v)));
        }
    }

    fn check_forest_state(board: &Board, state: &ForestState) {
        let active = board
            .dryness
            .iter()
            .enumerate()
            .map(|(v, &d)| d > 0 && !state.cuts[v])
            .collect();
        let mut rebuilt = split_tree_sets(board, active);
        let mut cached = Vec::new();
        for (id, set) in state.sets.iter().enumerate() {
            if let Some(set) = set {
                for &v in &set.cells {
                    assert_eq!(state.labels[v], id);
                }
                cached.push(set.cells.clone());
            }
        }
        for sets in [&mut rebuilt, &mut cached] {
            for set in sets.iter_mut() {
                set.sort_unstable();
            }
            sets.sort_unstable();
        }
        assert_eq!(rebuilt, cached);
        for (v, &d) in board.dryness.iter().enumerate() {
            assert_eq!(state.labels[v] == NONE, d == 0 || state.cuts[v]);
        }
        assert_eq!(
            state.cut_count,
            state.cuts.iter().filter(|&&cut| cut).count()
        );
        let total: f64 = state.sets.iter().flatten().map(|set| set.expected).sum();
        assert!((state.expected - total).abs() < 1e-9);
    }

    #[test]
    fn difference_update_merges_splits_and_preserves_unaffected_sets() {
        let board = Board::new(&padded(&[&[9, 9, 9], &[0, 0, 0], &[9, 0, 0]]));
        let deadline = Instant::now() + Duration::from_secs(30);
        let mut state = ForestState::new(&board, &[1], deadline).unwrap();
        check_forest_state(&board, &state);
        let unaffected_id = state.labels[6];
        let unaffected_expected = state.sets[unaffected_id].as_ref().unwrap().expected;
        assert_ne!(state.labels[0], state.labels[2]);

        let ProposalResult::Changed(merge) = state.propose(&board, &[1], false, deadline) else {
            panic!("結合の近傍を生成できませんでした");
        };
        // 棄却しても元の状態に変更がない。
        check_forest_state(&board, &state);
        assert!(state.cuts[1]);
        state.accept(merge);
        assert_eq!(state.labels[0], state.labels[2]);
        check_forest_state(&board, &state);

        let ProposalResult::Changed(split) = state.propose(&board, &[1], true, deadline) else {
            panic!("分割の近傍を生成できませんでした");
        };
        state.accept(split);
        assert_ne!(state.labels[0], state.labels[2]);
        assert_eq!(state.labels[6], unaffected_id);
        assert_eq!(
            state.sets[unaffected_id].as_ref().unwrap().expected,
            unaffected_expected
        );
        check_forest_state(&board, &state);
    }

    #[test]
    fn overwrite_handles_mixed_states_and_empty_sets() {
        let board = Board::new(&padded(&[&[9, 9], &[9, 9]]));
        let deadline = Instant::now() + Duration::from_secs(30);
        let mut state = ForestState::new(&board, &[0, 1, 2, 3], deadline).unwrap();
        assert_eq!(state.expected, 0.0);
        for (selected, on) in [
            (&[0, 1][..], false),
            (&[1, 2][..], false),
            (&[0, 1, 2, 3][..], true),
            (&[0, 1, 2, 3][..], false),
        ] {
            let ProposalResult::Changed(proposal) = state.propose(&board, selected, on, deadline)
            else {
                panic!("上書きの近傍を生成できませんでした");
            };
            state.accept(proposal);
            check_forest_state(&board, &state);
        }
        assert!(matches!(
            state.propose(&board, &[0, 1], false, deadline),
            ProposalResult::Unchanged
        ));
        assert_eq!(state.cut_count, 0);
    }

    #[test]
    fn edge_swap_is_atomic_and_counts_different_edge_lengths() {
        let board = Board::new(&padded(&[&[9, 9, 9], &[9, 9, 9], &[9, 9, 9]]));
        let layout = Layout {
            areas: Vec::new(),
            edges: vec![vec![1, 4], vec![6, 7, 8], Vec::new()],
            candidates: vec![1, 4, 6, 7, 8],
        };
        let deadline = Instant::now() + Duration::from_secs(30);
        let mut state = ForestState::new(&board, &layout.edges[0], deadline).unwrap();
        let original_cuts = state.cuts.clone();
        let original_expected = state.expected;
        let mut rng = 0;
        assert_eq!(
            choose_edge_swap(&layout, &state.cuts, &mut rng),
            Some((0, 1))
        );
        let ProposalResult::Changed(proposal) = state.propose_overwrites(
            &board,
            &[(&layout.edges[0], false), (&layout.edges[1], true)],
            deadline,
        ) else {
            panic!("辺の入れ替えを生成できませんでした");
        };
        // 受理前にはどちらの辺も変更されていない。
        assert_eq!(state.cuts, original_cuts);
        assert_eq!(state.expected, original_expected);
        check_forest_state(&board, &state);
        state.accept(proposal);
        assert!(layout.edges[0].iter().all(|&v| !state.cuts[v]));
        assert!(layout.edges[1].iter().all(|&v| state.cuts[v]));
        assert_eq!(state.cut_count, 3);
        assert_eq!(
            choose_edge_swap(&layout, &state.cuts, &mut rng),
            Some((1, 0))
        );
        check_forest_state(&board, &state);
    }

    #[test]
    fn edge_swap_requires_both_on_and_off_nonempty_edges() {
        let layout = Layout {
            areas: Vec::new(),
            edges: vec![Vec::new(), vec![0, 1], vec![2, 3]],
            candidates: vec![0, 1, 2, 3],
        };
        let mut rng = 0;
        assert_eq!(choose_edge_swap(&layout, &[true; 4], &mut rng), None);
        assert_eq!(choose_edge_swap(&layout, &[false; 4], &mut rng), None);
        assert_eq!(
            choose_edge_swap(&layout, &[true, true, false, false], &mut rng),
            Some((1, 2))
        );
    }

    #[test]
    fn compact_sampling_uses_whole_board_lightning_probability() {
        let board = Board::new(&padded(&[&[5, 5], &[0, 0]]));
        let estimate = estimate_tree_set(
            &board,
            &[0, 1],
            40_000,
            123,
            Instant::now() + Duration::from_secs(30),
        )
        .unwrap();
        assert!((estimate - 0.375).abs() < 0.01);
        assert!(estimate_tree_set(&board, &[0, 1], 100, 123, Instant::now()).is_none());
    }

    #[test]
    fn schedules_follow_specified_endpoints() {
        assert_eq!(
            annealing_parameters(0.0),
            (START_TEMPERATURE, START_PENALTY)
        );
        let (temperature, penalty) = annealing_parameters(1.0);
        assert!((temperature - END_TEMPERATURE).abs() < 1e-12);
        assert!((penalty - END_PENALTY).abs() < 1e-12);
        let (temperature, penalty) = annealing_parameters(0.5);
        assert!((temperature - (START_TEMPERATURE + END_TEMPERATURE) / 2.0).abs() < 1e-12);
        assert!((penalty - (START_PENALTY * END_PENALTY).sqrt()).abs() < 1e-12);
    }
}