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; // 定数定義 pub const EXPECTATION_SAMPLES: usize = 1000; const HIGH_RISK_PERCENTS: [usize; 7] = [40, 50, 60, 70, 80, 90, 100]; const NONE: usize = usize::MAX; #[derive(Debug, PartialEq)] pub struct BurnEstimate { pub expected_burned: f64, /// 外周パディングを除いたN×N配列。添字は0始まりの盤面座標。 /// 木のないマス・伐採済みのマス、および集計対象外のマスは0。 pub burn_probabilities: Vec>, } fn next_expectation_seed() -> u64 { // 起動ごとにランダムな初期値を用意し、呼び出しごとにシードを進める。 static STATE: OnceLock = 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) } /// 各乾燥連結成分の「全体の本数×対象の本数」の和を返す。 /// 対象マスのburn_weightsには、そのマスが属する乾燥連結成分のサイズを加算。 /// dryの外周はfalseとし、訪問したマスをfalseに変更する。 fn component_target_sum( dry: &mut [bool], targets: &[bool], width: usize, stack: &mut Vec, burn_weights: &mut [u64], ) -> 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.clear(); stack.push(start); let mut cursor = 0; while cursor < stack.len() { let v = stack[cursor]; cursor += 1; for next in [v - width, v + width, v - 1, v + 1] { if dry[next] { dry[next] = false; stack.push(next); } } } let size = stack.len() as u64; let mut target_size = 0_u64; for &v in stack.iter() { if targets[v] { target_size += 1; burn_weights[v] += size; } } sum += size * target_size; } } sum } /// 期待焼失本数と各マスの焼失確率をモンテカルロ推定する。 /// drynessはmainと同じ、外周を0で埋めた(N+2)×(N+2)の盤面。 /// sample_countには通常EXPECTATION_SAMPLESを渡す。呼び出すたびにシードを変更する。 /// 各マスは乾燥連結成分Cに属するサンプルで|C|/N²の確率で焼ける。 /// 焼失確率の合計が期待焼失本数に等しい。 /// 計算量はO(sample_count * N²)、追加メモリはO(N²)。 pub fn estimate_expected_burned(dryness: &[Vec], sample_count: usize) -> BurnEstimate { estimate_expected_burned_with_seed(dryness, None, sample_count, next_expectation_seed()) } /// 指定集合に属する木の期待焼失本数と各マスの焼失確率を推定する。 /// treesは外周を含めない0始まりの座標(i, j)の集合。重複は1本として扱う。 /// 木のないマスや伐採済みのマスを含めても、その寄与は0になる。 /// 集合外の木も延焼を伝え、落雷位置は盤面全体から選ばれる。 /// 戻り値の焼失確率は指定集合のマスのみ格納し、対象外は0とする。 /// 各サンプルでΣ(|C_k| * |C_k ∩ trees|) / N²を計算して平均する。 /// 盤面表現・サンプル数・シード更新はestimate_expected_burnedと同じ。 /// 現時点では全盤面をサンプリング・探索するため、計算量はO(sample_count * N²)。 pub fn estimate_expected_burned_for_trees( dryness: &[Vec], trees: &[(usize, usize)], sample_count: usize, ) -> BurnEstimate { estimate_expected_burned_with_seed(dryness, Some(trees), sample_count, next_expectation_seed()) } // 固定シードの内部処理は、テストの再現性を保つために分離する。 fn estimate_expected_burned_with_seed( dryness: &[Vec], trees: Option<&[(usize, usize)]>, sample_count: usize, seed: u64, ) -> BurnEstimate { 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 burn_weights = vec![0_u64; width * width]; 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, &mut burn_weights) as f64; } let denominator = sample_count as f64 * (n * n) as f64; let burn_probabilities = (1..=n) .map(|i| { (1..=n) .map(|j| burn_weights[i * width + j] as f64 / denominator) .collect() }) .collect(); BurnEstimate { expected_burned: total / denominator, burn_probabilities, } } struct Board { n: usize, dryness: Vec, neighbors: Vec<[usize; 4]>, } impl Board { fn new(dryness: &[Vec]) -> 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, } } } #[derive(Clone, Copy, Debug, PartialEq, Eq)] pub enum CutDirection { X, Y, XPlusY, XMinusY, } impl CutDirection { fn name(self) -> &'static str { match self { Self::X => "X=a", Self::Y => "Y=a", Self::XPlusY => "X+Y=a", Self::XMinusY => "X-Y=a", } } fn project(self, cell: usize, n: usize) -> i32 { let x = (cell % n) as i32; let y = (cell / n) as i32; match self { Self::X => x, Self::Y => y, Self::XPlusY => x + y, Self::XMinusY => x - y, } } } #[derive(Debug)] pub struct BalancedCut { pub direction: CutDirection, pub a: i32, pub cells: Vec, pub left_count: usize, pub right_count: usize, } /// 各方向でaを全探索し、左右の要素数の差が最小の切断を求める。 /// その4候補の中で、線上の木の本数が最小のものを返す。 /// 切断線上の木は左右のどちらにも含めず、両側に木が残る切断のみ扱う。 /// 同じ不均衡なら短い切断、同じ長さなら不均衡が小さい切断を優先する。 pub fn find_balanced_cut(cells: &[usize], n: usize) -> Option { let mut best: Option = None; for direction in [ CutDirection::X, CutDirection::Y, CutDirection::XPlusY, CutDirection::XMinusY, ] { let values: Vec = cells.iter().map(|&v| direction.project(v, n)).collect(); let Some(&minimum) = values.iter().min() else { return None; }; let maximum = *values.iter().max().unwrap(); let mut direction_best: Option = None; for a in minimum..=maximum { let left_count = values.iter().filter(|&&v| v < a).count(); let right_count = values.iter().filter(|&&v| v > a).count(); if left_count == 0 || right_count == 0 { continue; } let cut_cells: Vec = cells .iter() .copied() .zip(&values) .filter_map(|(v, &value)| (value == a).then_some(v)) .collect(); if cut_cells.is_empty() { continue; } let candidate = BalancedCut { direction, a, cells: cut_cells, left_count, right_count, }; let key = (left_count.abs_diff(right_count), candidate.cells.len()); if direction_best.as_ref().map_or(true, |old| { key < (old.left_count.abs_diff(old.right_count), old.cells.len()) }) { direction_best = Some(candidate); } } if let Some(candidate) = direction_best { let key = ( candidate.cells.len(), candidate.left_count.abs_diff(candidate.right_count), ); if best.as_ref().map_or(true, |old| { key < (old.cells.len(), old.left_count.abs_diff(old.right_count)) }) { best = Some(candidate); } } } best } // activeで指定されたマスだけの4近傍連結成分を取り出す。 // 抽出対象外の木を経由する実際の延焼は、別途行う焼失確率の計算に含まれる。 fn split_tree_sets(board: &Board, mut active: Vec) -> Vec> { 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 } // 現在存在する木の上位high_risk_percent%だけを連結成分の計算対象にする。 // 件数は切り捨て、同確率なら行・列の順。無視した木も実盤面には残す。 fn high_risk_mask(board: &Board, estimate: &BurnEstimate, high_risk_percent: usize) -> Vec { assert!(high_risk_percent <= 100); let mut trees: Vec = (0..board.n * board.n) .filter(|&v| board.dryness[v] > 0) .collect(); trees.sort_unstable_by(|&a, &b| { estimate.burn_probabilities[b / board.n][b % board.n] .total_cmp(&estimate.burn_probabilities[a / board.n][a % board.n]) .then_with(|| a.cmp(&b)) }); let keep = trees.len() * high_risk_percent / 100; let mut active = vec![false; board.n * board.n]; for &v in &trees[..keep] { active[v] = true; } active } fn highest_expected_component( board: &Board, active: Vec, estimate: &BurnEstimate, ) -> Option<(Vec, f64)> { let mut best: Option<(Vec, f64)> = None; for cells in split_tree_sets(board, active) { let expected: f64 = cells .iter() .map(|&v| estimate.burn_probabilities[v / board.n][v % board.n]) .sum(); if best.as_ref().map_or(true, |(_, old)| expected > *old) { best = Some((cells, expected)); } } best } fn apply_cut(dryness: &mut [Vec], cut: &BalancedCut, cuts: &mut Vec, k: usize) -> bool { if cut.cells.len() > k - cuts.len() { return false; } let n = dryness.len() - 2; for &v in &cut.cells { assert!(dryness[v / n + 1][v % n + 1] > 0); dryness[v / n + 1][v % n + 1] = 0; cuts.push(v); } true } fn solve(dryness: &mut [Vec], k: usize, high_risk_percent: usize) -> Vec { let mut cuts = Vec::new(); while cuts.len() < k { // 無視している木も含め、伐採後の実盤面で毎回サンプリングし直す。 let estimate = estimate_expected_burned(dryness, EXPECTATION_SAMPLES); let board = Board::new(dryness); let active = high_risk_mask(&board, &estimate, high_risk_percent); let Some((cells, component_expected)) = highest_expected_component(&board, active, &estimate) else { break; }; let Some(cut) = find_balanced_cut(&cells, board.n) else { eprintln!( "{:.6}\t{high_risk_percent}\t{component_expected:.6}\t{}\tNA\tNA\tNA\tNA\t0\t{}\tunsplittable", estimate.expected_burned, cells.len(), cuts.len() ); break; }; let applied = apply_cut(dryness, &cut, &mut cuts, k); let status = if applied { "cut" } else { "budget_limit" }; eprintln!( "{:.6}\t{high_risk_percent}\t{component_expected:.6}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{status}", estimate.expected_burned, cells.len(), cut.direction.name(), cut.a, cut.left_count, cut.right_count, cut.cells.len(), cuts.len() ); if !applied { break; } } cuts } struct Solution { high_risk_percent: usize, cuts: Vec, expected_burned: f64, } fn solve_best(dryness: &[Vec], k: usize) -> Solution { eprintln!("expected_burned\thigh_risk_percent\tcomponent_expected\tcomponent_size\tdirection\ta\tleft\tright\tcut_count\ttotal_cuts\tstatus"); let mut solutions = Vec::with_capacity(HIGH_RISK_PERCENTS.len()); for high_risk_percent in HIGH_RISK_PERCENTS { // 全候補を同じ初期盤面から開始し、前の候補の伐採を引き継がない。 let mut board = dryness.to_vec(); let cuts = solve(&mut board, k, high_risk_percent); // 最後の切断後の盤面を評価して比較する。 let expected_burned = estimate_expected_burned(&board, EXPECTATION_SAMPLES).expected_burned; solutions.push(Solution { high_risk_percent, cuts, expected_burned, }); } let best_index = (0..solutions.len()) .min_by(|&a, &b| { solutions[a] .expected_burned .total_cmp(&solutions[b].expected_burned) .then_with(|| { solutions[a] .high_risk_percent .cmp(&solutions[b].high_risk_percent) }) }) .unwrap(); eprintln!("\nexpected_burned\thigh_risk_percent\tcuts\tselected"); for (index, solution) in solutions.iter().enumerate() { eprintln!( "{:.6}\t{}\t{}\t{}", solution.expected_burned, solution.high_risk_percent, solution.cuts.len(), index == best_index ); } solutions.swap_remove(best_index) } fn main() { 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 cuts = solve_best(&dryness, k).cuts; println!("{}", cuts.len()); for v in cuts { println!("{} {}", v / n, v % n); } } #[cfg(test)] mod tests { use super::*; fn padded(rows: &[&[u32]]) -> Vec> { 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() { let result = estimate_expected_burned(&padded(&[&[0]]), 10); assert_eq!(result.expected_burned, 0.0); assert_eq!(result.burn_probabilities, vec![vec![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; } let mut weights = vec![0; 25]; assert_eq!( component_target_sum(&mut dry, &[true; 25], 5, &mut Vec::new(), &mut weights), 5 ); assert_eq!(weights[6], 2); assert_eq!(weights[7], 2); assert_eq!(weights[13], 1); 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.expected_burned - 0.7).abs() < 0.01); assert_eq!(estimate.burn_probabilities[0][0], estimate.expected_burned); 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.expected_burned - 0.375).abs() < 0.01); for probability in &estimate.burn_probabilities[0] { assert!((probability - 0.1875).abs() < 0.01); } assert_eq!(estimate.burn_probabilities[1], vec![0.0, 0.0]); let sum: f64 = estimate.burn_probabilities.iter().flatten().sum(); assert!((sum - estimate.expected_burned).abs() < 1e-12); } #[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.expected_burned - 0.1875).abs() < 0.01); assert_eq!(estimate.burn_probabilities[0][0], estimate.expected_burned); assert_eq!(estimate.burn_probabilities[0][1], 0.0); } #[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).expected_burned; 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).expected_burned, 0.0 ); assert_eq!( estimate_expected_burned_for_trees(&board, &[(0, 1)], 10).expected_burned, 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 balanced_cut_chooses_shorter_axis_for_rectangle() { let n = 5; let cells: Vec = (0..3) .flat_map(|i| (0..5).map(move |j| i * n + j)) .collect(); let cut = find_balanced_cut(&cells, n).unwrap(); assert_eq!(cut.direction, CutDirection::X); assert_eq!(cut.a, 2); assert_eq!(cut.cells.len(), 3); assert_eq!((cut.left_count, cut.right_count), (6, 6)); assert_eq!( cut.left_count + cut.right_count + cut.cells.len(), cells.len() ); } #[test] fn balanced_cut_uses_diagonal_and_excludes_cut_cells_from_halves() { let n = 7; let cells: Vec = (0_usize..n * n) .filter(|&v| (v / n).abs_diff(3) + (v % n).abs_diff(3) <= 2) .collect(); let cut = find_balanced_cut(&cells, n).unwrap(); assert_eq!(cut.direction, CutDirection::XPlusY); assert_eq!(cut.a, 6); assert_eq!(cut.cells.len(), 3); assert_eq!((cut.left_count, cut.right_count), (5, 5)); let cut = find_balanced_cut(&[0, 1, 2, 3], 2).unwrap(); assert_eq!(cut.cells.len(), 2); assert_eq!((cut.left_count, cut.right_count), (1, 1)); assert_eq!(CutDirection::XMinusY.project(3 * 5 + 1, 5), -2); assert!(find_balanced_cut(&[0, 1], 2).is_none()); } #[test] fn high_risk_selection_excludes_removed_trees_and_reorders_each_time() { let mut dryness = vec![vec![0; 7]; 7]; for row in &mut dryness[1..6] { row[1..6].fill(1); } let mut estimate = BurnEstimate { expected_burned: 0.0, burn_probabilities: vec![vec![0.5; 5]; 5], }; let mask = high_risk_mask(&Board::new(&dryness), &estimate, 20); assert_eq!(mask.iter().filter(|&&v| v).count(), 5); assert!(mask[..5].iter().all(|&v| v)); dryness[1][1] = 0; estimate.burn_probabilities[4].fill(0.9); let mask = high_risk_mask(&Board::new(&dryness), &estimate, 20); assert_eq!(mask.iter().filter(|&&v| v).count(), 4); assert!(!mask[0]); assert!(mask[20..24].iter().all(|&v| v)); } #[test] fn component_selection_uses_expected_burn_count_not_size() { let board = Board::new(&padded(&[&[9, 9, 0], &[0, 0, 0], &[0, 0, 9]])); let estimate = BurnEstimate { expected_burned: 0.7, burn_probabilities: vec![vec![0.1, 0.1, 0.0], vec![0.0; 3], vec![0.0, 0.0, 0.5]], }; let active = board.dryness.iter().map(|&d| d > 0).collect(); let (cells, expected) = highest_expected_component(&board, active, &estimate).unwrap(); assert_eq!(cells, vec![8]); assert_eq!(expected, 0.5); } #[test] fn cut_over_budget_is_not_partially_applied() { let mut board = padded(&[&[9, 9], &[9, 9]]); let original = board.clone(); let cut = find_balanced_cut(&[0, 1, 2, 3], 2).unwrap(); let mut cuts = Vec::new(); assert!(!apply_cut(&mut board, &cut, &mut cuts, 1)); assert_eq!(board, original); assert!(cuts.is_empty()); assert!(apply_cut(&mut board, &cut, &mut cuts, 2)); assert_eq!(cuts.len(), 2); for v in cuts { assert_eq!(board[v / 2 + 1][v % 2 + 1], 0); } } }