結果
| 問題 | No.5025 Forest Fire |
| コンテスト | |
| ユーザー |
👑  terry_u16
|
| 提出日時 | 2026-10-25 09:50:41 |
| 言語 | Rust (1.99.0 + proconio + num + itertools + ACL) |
| 結果 |
AC
不安定
|
| 実行時間 | 1,975 ms / 2,000 ms |
| + 950µs | |
| コード長 | 127,717 bytes |
| 記録 | |
| コンパイル時間 | 9,627 ms |
| コンパイル使用メモリ | 226,308 KB |
| 実行使用メモリ | 47,492 KB |
| スコア | 200,391,893 |
| 最終ジャッジ日時 | 2026-10-25 09:53:31 |
| 合計ジャッジ時間 | 169,367 ms |
|
ジャッジサーバーID （参考情報） |
judge5_0 / <nil> |
| 純コード判定しない問題か言語 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| other | AC * 50 |
ソースコード
pub mod library {
//! ヒューリスティック問題向けの共通ライブラリ。
pub mod bitboard {
use crate::library::grid::Coord;
/// 一辺64マス以下の正方形グリッドを、1行1ワードで表すビット集合。
#[derive(Debug, Clone, PartialEq, Eq)]
pub struct BitBoard<const N: usize> {
rows: [u64; N],
}
impl<const N: usize> Default for BitBoard<N> {
fn default() -> Self {
assert!(N > 0 && N <= 64);
Self { rows: [0; N] }
}
}
impl<const N: usize> BitBoard<N> {
pub fn insert(&mut self, p: Coord) {
debug_assert!(p.row() < N && p.col() < N);
self.rows[p.row()] |= 1_u64 << p.col();
}
pub fn contains(&self, p: Coord) -> bool {
debug_assert!(p.row() < N && p.col() < N);
self.rows[p.row()] & (1_u64 << p.col()) != 0
}
/// この集合の全マスがotherにも含まれるかを返す。
pub fn is_subset(&self, other: &Self) -> bool {
self.rows
.iter()
.zip(other.rows.iter())
.all(|(&a, &b)| a & !b == 0)
}
/// この集合内だけを4近傍移動し、startから距離radius以下のマスを返す。
/// 配列はスタック上に置き、毎段のヒープ確保・頂点キュー操作を避ける。
pub fn within_distance(&self, start: Coord, radius: usize) -> Self {
let mut reached = Self::default();
if !self.contains(start) {
return reached;
}
reached.insert(start);
let mut frontier = reached.rows;
let mut first = start.row();
let mut last = start.row();
for _ in 0..radius {
first = first.saturating_sub(1);
last = (last + 1).min(N - 1);
let mut next = [0_u64; N];
let mut any = 0;
// radius段で到達し得る行だけを走査する。
for row in first..=last {
let mut adjacent = (frontier[row] << 1) | (frontier[row] >> 1);
if row > 0 {
adjacent |= frontier[row - 1];
}
if row + 1 < N {
adjacent |= frontier[row + 1];
}
// 成分外・行の右端より外・既訪問のビットを落とす。
let bits = adjacent & self.rows[row] & !reached.rows[row];
next[row] = bits;
reached.rows[row] |= bits;
any |= bits;
}
if any == 0 {
break;
}
frontier = next;
}
reached
}
}
/// 行内の連続区間を一括で広げる、連結成分サイズ二乗和の作業領域。
pub struct ComponentWorkspace<const N: usize> {
queue: Vec<usize>,
}
impl<const N: usize> Default for ComponentWorkspace<N> {
fn default() -> Self {
Self {
queue: Vec::with_capacity(N * N),
}
}
}
impl<const N: usize> ComponentWorkspace<N> {
pub fn square_sum(&mut self, board: &BitBoard<N>, removed: &BitBoard<N>) -> u64 {
#[cfg(target_arch = "x86_64")]
if std::arch::is_x86_feature_detected!("popcnt") {
// 機能検出後だけ呼び、POPCNTのないCPUでは通常実装にフォールバックする。
return unsafe { self.square_sum_popcnt(board, removed) };
}
self.square_sum_impl(board, removed)
}
#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "popcnt")]
unsafe fn square_sum_popcnt(
&mut self,
board: &BitBoard<N>,
removed: &BitBoard<N>,
) -> u64 {
self.square_sum_impl(board, removed)
}
#[inline(always)]
fn square_sum_impl(&mut self, board: &BitBoard<N>, removed: &BitBoard<N>) -> u64 {
let mut active =
std::array::from_fn::<_, N, _>(|row| board.rows[row] & !removed.rows[row]);
let mut pending = [0_u64; N];
let mut total = 0;
for row in 0..N {
while active[row] != 0 {
pending[row] = active[row] & active[row].wrapping_neg();
self.queue.clear();
self.queue.push(row);
let mut head = 0;
let mut size = 0_u64;
while head < self.queue.len() {
let r = self.queue[head];
head += 1;
let mut reached = pending[r];
pending[r] = 0;
loop {
let expanded =
reached | (((reached << 1) | (reached >> 1)) & active[r]);
if expanded == reached {
break;
}
reached = expanded;
}
active[r] &= !reached;
size += u64::from(reached.count_ones());
for next in [r.wrapping_sub(1), r + 1] {
if next >= N {
continue;
}
let bits = reached & active[next];
if bits == 0 {
continue;
}
if pending[next] == 0 {
self.queue.push(next);
}
pending[next] |= bits;
}
}
total += size * size;
}
}
total
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::library::{
grid::{ADJACENTS, Map2d},
random::Random,
};
#[test]
fn subset_checks_empty_and_highest_bit() {
let empty = BitBoard::<64>::default();
let mut one = empty.clone();
one.insert(Coord::new(63, 63));
let mut two = one.clone();
two.insert(Coord::new(0, 0));
assert!(empty.is_subset(&one));
assert!(one.is_subset(&two));
assert!(two.is_subset(&two));
assert!(!one.is_subset(&empty));
assert!(!two.is_subset(&one));
}
#[test]
fn component_sums_include_bit63_and_workspace_reuse() {
let mut board = BitBoard::<64>::default();
for col in 0..64 {
board.insert(Coord::new(0, col));
board.insert(Coord::new(63, col));
}
let mut removed = BitBoard::default();
let mut workspace = ComponentWorkspace::default();
assert_eq!(workspace.square_sum(&board, &removed), 2 * 64 * 64);
removed.insert(Coord::new(0, 31));
assert_eq!(
workspace.square_sum(&board, &removed),
31 * 31 + 32 * 32 + 64 * 64
);
removed.insert(Coord::new(63, 63));
assert_eq!(
workspace.square_sum(&board, &removed),
31 * 31 + 32 * 32 + 63 * 63
);
assert_eq!(
workspace.square_sum_impl(&board, &removed),
31 * 31 + 32 * 32 + 63 * 63
);
assert_eq!(workspace.square_sum(&BitBoard::default(), &removed), 0);
}
fn queue_bfs<const N: usize>(
board: &BitBoard<N>,
start: Coord,
radius: usize,
) -> BitBoard<N> {
let mut reached = BitBoard::default();
if !board.contains(start) {
return reached;
}
let mut distances = Map2d::from_fn(|_| usize::MAX, N);
let mut queue = vec![start];
distances[start] = 0;
reached.insert(start);
let mut head = 0;
while head < queue.len() {
let p = queue[head];
head += 1;
if distances[p] == radius {
continue;
}
for diff in ADJACENTS {
let next = p + diff;
if next.in_map(N) && board.contains(next) && distances[next] == usize::MAX {
distances[next] = distances[p] + 1;
reached.insert(next);
queue.push(next);
}
}
}
reached
}
#[test]
fn detour_uses_path_distance_not_manhattan_distance() {
let mut board = BitBoard::<5>::default();
for row in 0..5 {
board.insert(Coord::new(row, 0));
board.insert(Coord::new(row, 2));
}
board.insert(Coord::new(4, 1));
let start = Coord::new(0, 0);
let goal = Coord::new(0, 2);
assert_eq!(start.dist(&goal), 2);
assert!(!board.within_distance(start, 9).contains(goal));
assert!(board.within_distance(start, 10).contains(goal));
}
#[test]
fn random_and_full_boards_match_queue_bfs() {
let mut rng = Random::new(123);
for density in [0, 4, 7, 10] {
for _ in 0..8 {
let mut board = BitBoard::<50>::default();
for row in 0..50 {
for col in 0..50 {
if rng.index(10) < density {
board.insert(Coord::new(row, col));
}
}
}
for start in [Coord::new(0, 0), Coord::new(49, 49), Coord::new(25, 25)] {
for radius in [0, 1, 4, 12, 60] {
assert_eq!(
board.within_distance(start, radius),
queue_bfs(&board, start, radius)
);
}
}
}
}
}
#[test]
fn edges_do_not_wrap_between_rows_or_columns() {
let mut board = BitBoard::<64>::default();
for p in [
Coord::new(0, 0),
Coord::new(0, 63),
Coord::new(1, 0),
Coord::new(63, 63),
] {
board.insert(p);
}
let start = Coord::new(0, 63);
let reached = board.within_distance(start, 100);
let mut expected = BitBoard::default();
expected.insert(start);
assert_eq!(reached, expected);
assert_eq!(reached, queue_bfs(&board, start, 100));
}
}
}
pub mod diagnostics {
use std::time::{Duration, Instant};
/// プログラムの経過時間を計測する。
pub struct Timer {
start: Instant,
}
impl Default for Timer {
fn default() -> Self {
Self {
start: Instant::now(),
}
}
}
impl Timer {
pub fn elapsed(&self) -> Duration {
self.start.elapsed()
}
}
}
pub mod grid {
use std::{
fmt::Display,
ops::{
Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Neg, Sub,
SubAssign,
},
};
#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash, Default)]
pub struct Coord {
row: u8,
col: u8,
}
impl Coord {
pub const fn new(row: usize, col: usize) -> Self {
Self {
row: row as u8,
col: col as u8,
}
}
pub const fn row(&self) -> usize {
self.row as usize
}
pub const fn col(&self) -> usize {
self.col as usize
}
pub fn in_map(&self, size: usize) -> bool {
self.row < size as u8 && self.col < size as u8
}
pub const fn to_index(&self, size: usize) -> CoordIndex {
CoordIndex(self.row as usize * size + self.col as usize)
}
pub const fn dist(&self, other: &Self) -> usize {
Self::dist_1d(self.row, other.row) + Self::dist_1d(self.col, other.col)
}
const fn dist_1d(x0: u8, x1: u8) -> usize {
x0.abs_diff(x1) as usize
}
}
impl Display for Coord {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
write!(f, "({}, {})", self.row, self.col)
}
}
#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash, Default)]
pub struct CoordIndex(pub usize);
impl CoordIndex {
pub const fn to_coord(&self, size: usize) -> Coord {
Coord::new(self.0 / size, self.0 % size)
}
}
#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash, Default)]
pub struct CoordDiff {
dr: i8,
dc: i8,
}
impl CoordDiff {
pub const fn new(dr: i32, dc: i32) -> Self {
Self {
dr: dr as i8,
dc: dc as i8,
}
}
pub const fn invert(&self) -> Self {
Self {
dr: -self.dr,
dc: -self.dc,
}
}
pub const fn dr(&self) -> i32 {
self.dr as i32
}
pub const fn dc(&self) -> i32 {
self.dc as i32
}
}
impl Display for CoordDiff {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
write!(f, "({}, {})", self.dr, self.dc)
}
}
impl Add<CoordDiff> for Coord {
type Output = Coord;
fn add(self, rhs: CoordDiff) -> Self::Output {
Coord {
row: self.row.wrapping_add(rhs.dr as u8),
col: self.col.wrapping_add(rhs.dc as u8),
}
}
}
impl AddAssign<CoordDiff> for Coord {
fn add_assign(&mut self, rhs: CoordDiff) {
self.row = self.row.wrapping_add(rhs.dr as u8);
self.col = self.col.wrapping_add(rhs.dc as u8);
}
}
impl Add<CoordDiff> for CoordDiff {
type Output = CoordDiff;
fn add(self, rhs: CoordDiff) -> Self::Output {
CoordDiff {
dr: self.dr + rhs.dr,
dc: self.dc + rhs.dc,
}
}
}
impl Add<Coord> for CoordDiff {
type Output = Coord;
fn add(self, rhs: Coord) -> Self::Output {
rhs + self
}
}
impl Sub<CoordDiff> for Coord {
type Output = Coord;
fn sub(self, rhs: CoordDiff) -> Self::Output {
Coord {
row: self.row.wrapping_sub(rhs.dr as u8),
col: self.col.wrapping_sub(rhs.dc as u8),
}
}
}
impl SubAssign<CoordDiff> for Coord {
fn sub_assign(&mut self, rhs: CoordDiff) {
*self = *self - rhs;
}
}
// Differences outside the i8 range wrap, as in CoordDiff::new.
impl Sub<Coord> for Coord {
type Output = CoordDiff;
fn sub(self, rhs: Coord) -> Self::Output {
CoordDiff::new(
self.row as i32 - rhs.row as i32,
self.col as i32 - rhs.col as i32,
)
}
}
impl Sub<CoordDiff> for CoordDiff {
type Output = CoordDiff;
fn sub(self, rhs: CoordDiff) -> Self::Output {
CoordDiff {
dr: self.dr - rhs.dr,
dc: self.dc - rhs.dc,
}
}
}
impl SubAssign<CoordDiff> for CoordDiff {
fn sub_assign(&mut self, rhs: CoordDiff) {
*self = *self - rhs;
}
}
// Scalar operations use i32 arithmetic, then narrow as in CoordDiff::new.
impl Mul<i32> for CoordDiff {
type Output = CoordDiff;
fn mul(self, rhs: i32) -> Self::Output {
CoordDiff::new(self.dr() * rhs, self.dc() * rhs)
}
}
impl Mul<CoordDiff> for i32 {
type Output = CoordDiff;
fn mul(self, rhs: CoordDiff) -> Self::Output {
rhs * self
}
}
impl MulAssign<i32> for CoordDiff {
fn mul_assign(&mut self, rhs: i32) {
*self = *self * rhs;
}
}
impl Div<i32> for CoordDiff {
type Output = CoordDiff;
fn div(self, rhs: i32) -> Self::Output {
CoordDiff::new(self.dr() / rhs, self.dc() / rhs)
}
}
impl DivAssign<i32> for CoordDiff {
fn div_assign(&mut self, rhs: i32) {
*self = *self / rhs;
}
}
impl AddAssign<CoordDiff> for CoordDiff {
fn add_assign(&mut self, rhs: CoordDiff) {
self.dr += rhs.dr;
self.dc += rhs.dc;
}
}
impl Neg for CoordDiff {
type Output = Self;
fn neg(self) -> Self::Output {
Self {
dr: -self.dr,
dc: -self.dc,
}
}
}
pub const ADJACENTS: [CoordDiff; 4] = [
CoordDiff::new(-1, 0),
CoordDiff::new(0, 1),
CoordDiff::new(1, 0),
CoordDiff::new(0, -1),
];
pub const U: usize = 0;
pub const R: usize = 1;
pub const D: usize = 2;
pub const L: usize = 3;
pub const DIRECTIONS: [char; 4] = ['U', 'R', 'D', 'L'];
#[derive(Debug, Clone)]
pub struct Map2d<T> {
size: usize,
map: Vec<T>,
}
impl<T> Map2d<T> {
pub fn new(map: Vec<T>, size: usize) -> Self {
debug_assert!(size * size == map.len());
Self { size, map }
}
pub fn from_fn(mut f: impl FnMut(Coord) -> T, size: usize) -> Self {
let mut map = Vec::with_capacity(size * size);
for row in 0..size {
for col in 0..size {
map.push(f(Coord::new(row, col)));
}
}
Self { size, map }
}
pub fn iter(&self) -> impl Iterator<Item = &T> {
self.map.iter()
}
}
impl<T: Default + Clone> Map2d<T> {
pub fn with_default(size: usize) -> Self {
let map = vec![T::default(); size * size];
Self { size, map }
}
}
impl<T> Index<Coord> for Map2d<T> {
type Output = T;
#[inline]
fn index(&self, coordinate: Coord) -> &Self::Output {
&self.map[coordinate.to_index(self.size).0]
}
}
impl<T> IndexMut<Coord> for Map2d<T> {
#[inline]
fn index_mut(&mut self, coordinate: Coord) -> &mut Self::Output {
&mut self.map[coordinate.to_index(self.size).0]
}
}
impl<T> Index<&Coord> for Map2d<T> {
type Output = T;
#[inline]
fn index(&self, coordinate: &Coord) -> &Self::Output {
&self.map[coordinate.to_index(self.size).0]
}
}
impl<T> IndexMut<&Coord> for Map2d<T> {
#[inline]
fn index_mut(&mut self, coordinate: &Coord) -> &mut Self::Output {
&mut self.map[coordinate.to_index(self.size).0]
}
}
impl<T> Index<usize> for Map2d<T> {
type Output = [T];
#[inline]
fn index(&self, row: usize) -> &Self::Output {
let begin = row * self.size;
let end = begin + self.size;
&self.map[begin..end]
}
}
impl<T> IndexMut<usize> for Map2d<T> {
#[inline]
fn index_mut(&mut self, row: usize) -> &mut Self::Output {
let begin = row * self.size;
let end = begin + self.size;
&mut self.map[begin..end]
}
}
impl<T> Index<CoordIndex> for Map2d<T> {
type Output = T;
#[inline]
fn index(&self, index: CoordIndex) -> &Self::Output {
&self.map[index.0]
}
}
impl<T> IndexMut<CoordIndex> for Map2d<T> {
#[inline]
fn index_mut(&mut self, index: CoordIndex) -> &mut Self::Output {
&mut self.map[index.0]
}
}
#[cfg(test)]
mod test {
use super::{Coord, CoordDiff, Map2d};
#[test]
fn coord_add() {
let c = Coord::new(2, 4);
let d = CoordDiff::new(-3, 5);
let actual = c + d;
let expected = Coord::new(!0, 9);
assert_eq!(expected, actual);
}
#[test]
fn coord_add_assign() {
let mut c = Coord::new(2, 4);
let d = CoordDiff::new(-3, 5);
c += d;
let expected = Coord::new(!0, 9);
assert_eq!(expected, c);
}
#[test]
fn coord_translation_round_trip() {
for row in 0..=255 {
let c = Coord::new(row, 255 - row);
for dr in -128..=127 {
let d = CoordDiff::new(dr, -1 - dr);
assert_eq!(c + d, d + c);
assert_eq!((c + d) - d, c);
assert_eq!((c + d) - c, d);
let mut translated = c;
translated -= d;
assert_eq!(translated + d, c);
}
}
}
#[test]
fn coord_difference() {
assert_eq!(Coord::new(2, 9) - Coord::new(5, 4), CoordDiff::new(-3, 5));
assert_eq!(
Coord::new(255, 0) - Coord::new(0, 255),
CoordDiff::new(-1, 1)
);
}
#[test]
fn coord_diff_arithmetic() {
let a = CoordDiff::new(-3, 5);
let b = CoordDiff::new(7, -2);
assert_eq!(a + b, CoordDiff::new(4, 3));
assert_eq!(a - b, CoordDiff::new(-10, 7));
assert_eq!(-a, a.invert());
assert_eq!(a + (-a), CoordDiff::default());
let mut d = a;
d += b;
assert_eq!(d, a + b);
d -= b;
assert_eq!(d, a);
}
#[test]
fn coord_diff_scalar_arithmetic() {
let d = CoordDiff::new(-3, 5);
assert_eq!(d * -2, CoordDiff::new(6, -10));
assert_eq!(-2 * d, d * -2);
assert_eq!(d * 0, CoordDiff::default());
assert_eq!(d / 2, CoordDiff::new(-1, 2));
assert_eq!(d / -2, CoordDiff::new(1, -2));
assert_eq!(CoordDiff::new(64, -128) * 2, CoordDiff::new(-128, 0));
assert_eq!(CoordDiff::new(-128, 0) / -1, CoordDiff::new(-128, 0));
assert_eq!(d / 256, CoordDiff::default());
let mut scaled = d;
scaled *= -2;
assert_eq!(scaled, CoordDiff::new(6, -10));
scaled /= -2;
assert_eq!(scaled, d);
}
#[test]
#[should_panic]
fn coord_diff_divide_by_zero() {
let _ = CoordDiff::new(1, 2) / 0;
}
#[test]
fn map_new() {
let map = Map2d::new(vec![0, 1, 2, 3], 2);
let actual = map[Coord::new(1, 0)];
let expected = 2;
assert_eq!(expected, actual);
}
#[test]
fn map_default() {
let map = Map2d::with_default(2);
let actual = map[Coord::new(1, 0)];
let expected = 0;
assert_eq!(expected, actual);
}
}
}
pub mod maxflow {
//! 辺を追加してから解く、CSR隣接配列による整数最大流。
#[derive(Clone, Copy)]
struct Edge {
to: u16,
capacity: u32,
}
/// 頂点数は65,535以下、各辺の容量はu32以内。flow後の辺追加は行わない。
pub struct MaxFlow {
edges: Vec<Edge>,
degree: Vec<usize>,
start: Vec<usize>,
adjacency: Vec<u32>,
level: Vec<i32>,
cursor: Vec<usize>,
queue: Vec<usize>,
}
impl MaxFlow {
pub fn new(n: usize) -> Self {
assert!(n > 0 && n <= u16::MAX as usize);
Self {
edges: Vec::with_capacity(n * 8),
degree: vec![0; n],
start: Vec::new(),
adjacency: Vec::new(),
level: vec![-1; n],
cursor: vec![0; n],
queue: Vec::with_capacity(n),
}
}
pub fn add_edge(&mut self, from: usize, to: usize, capacity: i64) -> usize {
assert!(self.start.is_empty());
assert!(from < self.degree.len() && to < self.degree.len());
let capacity = u32::try_from(capacity).unwrap();
let id = self.edges.len() / 2;
self.edges.push(Edge {
to: to as u16,
capacity,
});
self.edges.push(Edge {
to: from as u16,
capacity: 0,
});
self.degree[from] += 1;
self.degree[to] += 1;
id
}
fn prepare(&mut self) {
if !self.start.is_empty() {
return;
}
let n = self.degree.len();
self.start = vec![0; n + 1];
for v in 0..n {
self.start[v + 1] = self.start[v] + self.degree[v];
}
self.adjacency = vec![0; self.edges.len()];
self.cursor.copy_from_slice(&self.start[..n]);
for i in 0..self.edges.len() {
let from = self.edges[i ^ 1].to as usize;
self.adjacency[self.cursor[from]] = i as u32;
self.cursor[from] += 1;
}
}
fn bfs(&mut self, source: usize, sink: usize) {
self.level.fill(-1);
self.level[source] = 0;
self.queue.clear();
self.queue.push(source);
let mut head = 0;
while head < self.queue.len() {
let v = self.queue[head];
head += 1;
for &id in &self.adjacency[self.start[v]..self.start[v + 1]] {
let edge = self.edges[id as usize];
let to = edge.to as usize;
if edge.capacity == 0 || self.level[to] >= 0 {
continue;
}
self.level[to] = self.level[v] + 1;
if to == sink {
return;
}
self.queue.push(to);
}
}
}
fn dfs(&mut self, v: usize, source: usize, limit: u32) -> u32 {
if v == source {
return limit;
}
let mut sent = 0;
let level = self.level[v];
while self.cursor[v] < self.start[v + 1] {
let id = self.adjacency[self.cursor[v]] as usize;
let to = self.edges[id].to as usize;
let capacity = self.edges[id ^ 1].capacity;
if level > self.level[to] && capacity > 0 {
let amount = self.dfs(to, source, (limit - sent).min(capacity));
self.edges[id].capacity += amount;
self.edges[id ^ 1].capacity -= amount;
sent += amount;
if sent == limit {
return sent;
}
}
self.cursor[v] += 1;
}
sent
}
pub fn flow(&mut self, source: usize, sink: usize) -> u64 {
assert_ne!(source, sink);
assert!(source < self.degree.len() && sink < self.degree.len());
self.prepare();
let mut result = 0;
loop {
self.bfs(source, sink);
if self.level[sink] < 0 {
return result;
}
self.cursor
.copy_from_slice(&self.start[..self.degree.len()]);
loop {
let sent = self.dfs(sink, source, u32::MAX);
if sent == 0 {
break;
}
result += u64::from(sent);
}
}
}
pub fn min_cut(&self, source: usize) -> Vec<bool> {
assert!(!self.start.is_empty());
let mut visited = vec![false; self.degree.len()];
let mut queue = vec![source];
visited[source] = true;
let mut head = 0;
while head < queue.len() {
let v = queue[head];
head += 1;
for &id in &self.adjacency[self.start[v]..self.start[v + 1]] {
let edge = self.edges[id as usize];
let to = edge.to as usize;
if edge.capacity > 0 && !visited[to] {
visited[to] = true;
queue.push(to);
}
}
}
visited
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::library::random::Random;
#[test]
fn random_networks_match_acl_flow_and_cut() {
let mut rng = Random::new(2026);
for _ in 0..1000 {
let n = 2 + rng.index(30);
let mut graph = MaxFlow::new(n);
let mut reference = ac_library::MfGraph::<i64>::new(n);
for _ in 0..n * 10 {
let from = rng.index(n);
let to = rng.index(n);
let capacity = rng.index(1000) as i64;
graph.add_edge(from, to, capacity);
reference.add_edge(from, to, capacity);
}
assert_eq!(graph.flow(0, n - 1), reference.flow(0, n - 1) as u64);
assert_eq!(graph.min_cut(0), reference.min_cut(0));
assert_eq!(graph.flow(0, n - 1), 0);
}
}
#[test]
fn total_flow_can_exceed_single_edge_capacity_limit() {
let mut graph = MaxFlow::new(2);
graph.add_edge(0, 1, i64::from(u32::MAX));
graph.add_edge(0, 1, i64::from(u32::MAX));
assert_eq!(graph.flow(0, 1), u64::from(u32::MAX) * 2);
assert_eq!(graph.min_cut(0), [true, false]);
}
}
}
pub mod random {
/// 外部クレートに依存しないSplitMix64乱数生成器。
#[derive(Debug, Clone)]
pub struct Random {
state: u64,
}
impl Random {
pub fn new(seed: u64) -> Self {
Self { state: seed }
}
pub fn next_u64(&mut self) -> u64 {
self.state = self.state.wrapping_add(0x9e3779b97f4a7c15);
let mut z = self.state;
z = (z ^ (z >> 30)).wrapping_mul(0xbf58476d1ce4e5b9);
z = (z ^ (z >> 27)).wrapping_mul(0x94d049bb133111eb);
z ^ (z >> 31)
}
/// 棄却法により剰余の偏りを除き、0..upperから一様に選ぶ。
pub fn index(&mut self, upper: usize) -> usize {
assert!(upper > 0);
let bound = upper as u64;
let threshold = bound.wrapping_neg() % bound;
loop {
let value = self.next_u64();
if value >= threshold {
return (value % bound) as usize;
}
}
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn splitmix64_known_sequence() {
let mut rng = Random::new(0);
assert_eq!(rng.next_u64(), 0xe220a8397b1dcdaf);
assert_eq!(rng.next_u64(), 0x6e789e6aa1b965f4);
}
#[test]
fn index_stays_in_range() {
let mut rng = Random::new(42);
for upper in [1, 2, 10, 2500] {
for _ in 0..10000 {
assert!(rng.index(upper) < upper);
}
}
}
}
}
}
pub mod problem {
use std::fmt::{self, Display, Formatter};
use std::io::{self, BufRead};
use proconio::{input, marker::Bytes, source::once::OnceSource};
use crate::library::grid::{Coord, Map2d};
#[derive(Debug, Clone)]
pub struct Input {
/// 各マスの乾燥度。0は木のないマス、1〜9は木の乾燥度を表す。
pub d: Map2d<u8>,
}
impl Input {
/// 森の一辺のマス数。
pub const N: usize = 50;
/// 伐採できる木の本数の上限。
pub const K: usize = 100;
pub fn read() -> Self {
Self::read_from(io::stdin().lock())
}
fn read_from(reader: impl BufRead) -> Self {
let mut source = OnceSource::new(reader);
input! {
from &mut source,
n: usize,
k: usize,
}
assert_eq!(n, Self::N, "入力のNが定数と一致しません");
assert_eq!(k, Self::K, "入力のKが定数と一致しません");
let mut cells = Vec::with_capacity(Self::N * Self::N);
for _ in 0..Self::N {
input! { from &mut source, row: Bytes }
assert_eq!(row.len(), Self::N, "乾燥度の行の長さがNと一致しません");
for digit in row {
assert!(digit.is_ascii_digit(), "乾燥度は0〜9の数字で指定します");
cells.push(digit - b'0');
}
}
Self {
d: Map2d::new(cells, Self::N),
}
}
}
#[derive(Debug, Clone, Default)]
pub struct Output {
/// 伐採する木のマス（0始まりの行・列）。相異なる木のあるマスをK個以下指定する。
pub cuts: Vec<Coord>,
}
impl Display for Output {
fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
// 本数は座標列から求め、出力行数との不一致を防ぐ。
writeln!(f, "{}", self.cuts.len())?;
for coord in &self.cuts {
writeln!(f, "{} {}", coord.row(), coord.col())?;
}
Ok(())
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn reads_digits_in_row_major_order_including_leading_zero() {
let mut text = String::from("50 100\n");
for row in 0..50 {
for col in 0..50 {
text.push(char::from(b'0' + ((row + col) % 10) as u8));
}
text.push('\n');
}
let input = Input::read_from(text.as_bytes());
assert_eq!(input.d.iter().count(), 2500);
for row in 0..50 {
for col in 0..50 {
assert_eq!(input.d[Coord::new(row, col)], ((row + col) % 10) as u8);
}
}
}
#[test]
fn formats_sample_output() {
let output = Output {
cuts: vec![Coord::new(0, 32), Coord::new(1, 33), Coord::new(2, 34)],
};
assert_eq!(output.to_string(), "3\n0 32\n1 33\n2 34\n");
}
#[test]
fn formats_empty_output() {
assert_eq!(Output::default().to_string(), "0\n");
}
}
}
pub mod solver {
mod vertex_cut {
mod evaluation {
#[cfg(test)]
use crate::library::grid::CoordIndex;
use crate::library::{
bitboard::{BitBoard, ComponentWorkspace},
grid::{ADJACENTS, Coord, Map2d},
random::Random,
};
use crate::problem::Input;
const SAMPLES: usize = 128;
pub(super) struct Evaluator {
samples: Vec<Map2d<bool>>,
baseline: u64,
bit_samples: Vec<BitBoard<{ Input::N }>>,
pruned: Option<PrunedSamples>,
}
struct PrunedSamples {
allowed: BitBoard<{ Input::N }>,
samples: Vec<BitBoard<{ Input::N }>>,
constant: u64,
}
fn sample_seed(input: &Input) -> u64 {
// 配布ツールと同じ入力由来のseedで、先頭SAMPLES個のサンプルを固定する。
let mut hash = 1469598103934665603_u64;
for &d in input.d.iter() {
hash ^= u64::from(b'0' + d);
hash = hash.wrapping_mul(1099511628211);
}
hash ^ (Input::N as u64)
.wrapping_mul(1000003)
.wrapping_add((Input::K as u64).wrapping_mul(998244353))
}
// 周囲1マスを壁にして、BFS内の座標変換と盤面内判定を省く。
#[cfg(test)]
const PADDED_N: usize = Input::N + 2;
#[cfg(test)]
struct EvaluationWorkspace {
active: Map2d<bool>,
queue: Vec<usize>,
}
#[cfg(test)]
impl EvaluationWorkspace {
fn new() -> Self {
Self {
active: Map2d::with_default(PADDED_N),
queue: Vec::with_capacity(Input::N * Input::N),
}
}
fn component_square_sum(&mut self, dry: &Map2d<bool>, cuts: &[usize]) -> u64 {
for row in 0..Input::N {
self.active[row + 1][1..=Input::N].copy_from_slice(&dry[row]);
}
for &index in cuts {
self.active[CoordIndex(index)] = false;
}
let mut sum = 0;
for start in PADDED_N..(PADDED_N * (Input::N + 1)) {
if !self.active[CoordIndex(start)] {
continue;
}
self.active[CoordIndex(start)] = false;
self.queue.clear();
self.queue.push(start);
let mut head = 0;
while head < self.queue.len() {
let p = self.queue[head];
head += 1;
for next in [p - PADDED_N, p + 1, p + PADDED_N, p - 1] {
if self.active[CoordIndex(next)] {
self.active[CoordIndex(next)] = false;
self.queue.push(next);
}
}
}
sum += (self.queue.len() as u64).pow(2);
}
sum
}
}
#[cfg(test)]
fn padded_cuts(cut: &Map2d<bool>) -> Vec<usize> {
cut.iter()
.enumerate()
.filter_map(|(i, &removed)| {
removed.then_some((i / Input::N + 1) * PADDED_N + i % Input::N + 1)
})
.collect()
}
fn to_bitboard(map: &Map2d<bool>) -> BitBoard<{ Input::N }> {
let mut result = BitBoard::default();
for row in 0..Input::N {
for col in 0..Input::N {
if map[row][col] {
result.insert(crate::library::grid::Coord::new(row, col));
}
}
}
result
}
impl Evaluator {
pub(super) fn new(input: &Input) -> Self {
let mut this = Self::with_samples(input, SAMPLES);
this.prepare_pruned_samples(input);
this
}
pub(super) fn with_samples(input: &Input, count: usize) -> Self {
let mut rng = Random::new(sample_seed(input));
let thresholds: Vec<u64> = input
.d
.iter()
.map(|&d| u64::from(d) * (1_u64 << 32) / 10)
.collect();
let samples: Vec<_> = (0..count)
.map(|_| {
let mut cells = Vec::with_capacity(Input::N * Input::N);
for pair in thresholds.chunks(2) {
let bits = rng.next_u64();
cells.push((bits & 0xffff_ffff) < pair[0]);
if pair.len() == 2 {
cells.push((bits >> 32) < pair[1]);
}
}
Map2d::new(cells, Input::N)
})
.collect();
let bit_samples: Vec<_> = samples.iter().map(to_bitboard).collect();
let mut workspace = ComponentWorkspace::default();
let baseline = bit_samples
.iter()
.map(|sample| workspace.square_sum(sample, &BitBoard::default()))
.sum();
Self {
samples,
baseline,
bit_samples,
pruned: None,
}
}
fn prepare_pruned_samples(&mut self, input: &Input) {
let allowed = Map2d::from_fn(|p| input.d[p] >= super::D_THRESHOLD, Input::N);
let mut constant = 0;
let mut samples = Vec::with_capacity(self.samples.len());
let mut queue = Vec::with_capacity(Input::N * Input::N);
for sample in &self.samples {
let mut seen = Map2d::with_default(Input::N);
let mut retained = BitBoard::default();
for row in 0..Input::N {
for col in 0..Input::N {
let start = Coord::new(row, col);
if seen[start] || !sample[start] {
continue;
}
queue.clear();
queue.push(start);
seen[start] = true;
let mut head = 0;
let mut contains_allowed = false;
while head < queue.len() {
let p = queue[head];
head += 1;
contains_allowed |= allowed[p];
for diff in ADJACENTS {
let next = p + diff;
if next.in_map(Input::N) && !seen[next] && sample[next] {
seen[next] = true;
queue.push(next);
}
}
}
if contains_allowed {
for &p in &queue {
retained.insert(p);
}
} else {
constant += (queue.len() as u64).pow(2);
}
}
}
samples.push(retained);
}
self.pruned = Some(PrunedSamples {
allowed: to_bitboard(&allowed),
samples,
constant,
});
}
pub(super) fn samples(&self) -> &[Map2d<bool>] {
&self.samples
}
pub(super) fn score(&self, cut: &Map2d<bool>) -> i64 {
self.score_from_total(self.total(cut))
}
pub(super) fn total(&self, cut: &Map2d<bool>) -> u64 {
let cuts = to_bitboard(cut);
let mut workspace = ComponentWorkspace::default();
if let Some(pruned) = &self.pruned {
if cuts.is_subset(&pruned.allowed) {
return pruned.constant
+ pruned
.samples
.iter()
.map(|sample| workspace.square_sum(sample, &cuts))
.sum::<u64>();
}
}
self.bit_samples
.iter()
.map(|sample| workspace.square_sum(sample, &cuts))
.sum()
}
pub(super) fn score_from_total(&self, total: u64) -> i64 {
let denominator = (self.samples.len() * Input::N * Input::N) as f64;
let e0 = self.baseline as f64 / denominator;
let e = total as f64 / denominator;
// ゼロ近傍の処理・丸めを配布ツールと揃える。
(1e6 * (1.0 + (e0 / e.max(1e-12)).log2())).round() as i64
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::library::grid::{ADJACENTS, Coord};
#[test]
fn pruning_constant_components_preserves_scores_and_falls_back() {
let mut rng = Random::new(123);
let input = Input {
d: Map2d::from_fn(
|p| {
if p.row() < 15 && p.col() < 15 {
rng.index(10) as u8
} else {
0
}
},
Input::N,
),
};
let full = Evaluator::with_samples(&input, 8);
let mut pruned = Evaluator::with_samples(&input, 8);
pruned.prepare_pruned_samples(&input);
assert!(pruned.pruned.as_ref().unwrap().constant > 0);
for restricted in [true, false] {
for _ in 0..20 {
let cut = Map2d::from_fn(
|p| {
rng.index(5) == 0
&& (!restricted || input.d[p] >= super::super::D_THRESHOLD)
},
Input::N,
);
assert_eq!(full.total(&cut), pruned.total(&cut));
assert_eq!(full.score(&cut), pruned.score(&cut));
}
}
}
fn component_square_sum(dry: &Map2d<bool>, cut: &Map2d<bool>) -> u64 {
let mut seen = Map2d::with_default(Input::N);
let mut queue = Vec::with_capacity(Input::N * Input::N);
let mut sum = 0;
for row in 0..Input::N {
for col in 0..Input::N {
let start = Coord::new(row, col);
if seen[start] || !dry[start] || cut[start] {
continue;
}
seen[start] = true;
queue.clear();
queue.push(start);
let mut head = 0;
while head < queue.len() {
let p = queue[head];
head += 1;
for diff in ADJACENTS {
let next = p + diff;
if next.in_map(Input::N)
&& !seen[next]
&& dry[next]
&& !cut[next]
{
seen[next] = true;
queue.push(next);
}
}
}
sum += (queue.len() as u64).pow(2);
}
}
sum
}
#[test]
fn square_sum_and_cut_splitting() {
let mut dry = Map2d::with_default(Input::N);
for col in 0..3 {
dry[Coord::new(0, col)] = true;
}
dry[Coord::new(2, 2)] = true;
let mut cut = Map2d::with_default(Input::N);
let mut workspace = EvaluationWorkspace::new();
assert_eq!(workspace.component_square_sum(&dry, &padded_cuts(&cut)), 10);
cut[Coord::new(0, 1)] = true;
assert_eq!(workspace.component_square_sum(&dry, &padded_cuts(&cut)), 3);
}
#[test]
fn padded_bfs_matches_original_and_workspace_can_be_reused() {
let mut rng = Random::new(2026);
let mut workspace = EvaluationWorkspace::new();
let mut bit_workspace = ComponentWorkspace::default();
for density in [0, 1, 5, 9, 10] {
for _ in 0..20 {
let dry = Map2d::from_fn(|_| rng.index(10) < density, Input::N);
let cut = Map2d::from_fn(|_| rng.index(25) == 0, Input::N);
assert_eq!(
bit_workspace.square_sum(&to_bitboard(&dry), &to_bitboard(&cut)),
component_square_sum(&dry, &cut)
);
assert_eq!(
workspace.component_square_sum(&dry, &padded_cuts(&cut)),
component_square_sum(&dry, &cut)
);
}
}
}
#[test]
fn no_cut_score_is_one_million_and_samples_are_fixed() {
let input = Input {
d: Map2d::from_fn(|p| ((p.row() + p.col()) % 10) as u8, Input::N),
};
let first = Evaluator::new(&input);
let second = Evaluator::new(&input);
assert_eq!(first.score(&Map2d::with_default(Input::N)), 1_000_000);
for (a, b) in first.samples.iter().zip(&second.samples) {
assert!(a.iter().eq(b.iter()));
for row in 0..Input::N {
for col in 0..Input::N {
let p = Coord::new(row, col);
if input.d[p] == 0 {
assert!(!a[p]);
}
}
}
}
}
}
}
mod incremental {
//! 伐採の変更に関係する連結成分だけ、二乗和と全頂点削除利得を更新する。
use crate::library::grid::{Coord, CoordIndex, Map2d};
use crate::problem::Input;
pub(super) const N: usize = Input::N + 2;
pub(super) fn index(p: Coord) -> CoordIndex {
CoordIndex((p.row() + 1) * N + p.col() + 1)
}
// 一時領域は全サンプルで共用し、サンプルごとの確保・キャッシュ消費を避ける。
struct Workspace {
affected: Vec<u16>,
queue: Vec<u16>,
visited: Vec<(u16, u16, u32)>,
}
impl Workspace {
fn new() -> Self {
Self {
affected: Vec::with_capacity(Input::N * Input::N),
queue: Vec::with_capacity(Input::N * Input::N),
visited: Vec::with_capacity(Input::N * Input::N),
}
}
}
struct Sample {
dry: Map2d<bool>,
order: Map2d<u16>,
component: Map2d<u16>,
gain: Map2d<u32>,
// 成分IDは消滅時に回収する。全マス分の空Vecをサンプルごとに確保しない。
members: Vec<Vec<u16>>,
free_ids: Vec<usize>,
clock: u16,
total: u64,
}
impl Sample {
fn new(
dry: &Map2d<bool>,
cuts: &[Coord],
gains: &mut Map2d<u64>,
workspace: &mut Workspace,
) -> Self {
let mut this = Self {
dry: Map2d::with_default(N),
order: Map2d::from_fn(|_| u16::MAX, N),
component: Map2d::with_default(N),
gain: Map2d::with_default(N),
members: Vec::new(),
free_ids: Vec::new(),
clock: 0,
total: 0,
};
for row in 0..Input::N {
this.dry[row + 1][1..=Input::N].copy_from_slice(&dry[row]);
for col in 0..Input::N {
if dry[row][col] {
this.order[row + 1][col + 1] = 0;
}
}
}
for &p in cuts {
this.order[index(p)] = u16::MAX;
}
for root in N..N * (Input::N + 1) {
if this.order[CoordIndex(root)] == 0 {
this.rebuild_component(root, gains, workspace);
}
}
this
}
fn dfs(
&mut self,
v: usize,
parent: usize,
workspace: &mut Workspace,
) -> (u16, u32) {
self.clock += 1;
let order = self.clock;
self.order[CoordIndex(v)] = order;
let mut low = order;
let mut size = 1;
let mut separated_size = 1;
let mut separated_square = 0;
for next in [v - N, v + 1, v + N, v - 1] {
let next_order = self.order[CoordIndex(next)];
if next_order == 0 {
let (child_low, child_size) = self.dfs(next, v, workspace);
size += child_size;
low = low.min(child_low);
if child_low >= order {
separated_size += child_size;
separated_square += child_size * child_size;
}
} else if next != parent {
low = low.min(next_order);
}
}
workspace.visited.push((
v as u16,
separated_size as u16,
separated_size * separated_size + separated_square,
));
(low, size)
}
fn rebuild_component(
&mut self,
root: usize,
gains: &mut Map2d<u64>,
workspace: &mut Workspace,
) {
self.clock = 0;
workspace.visited.clear();
let (_, size) = self.dfs(root, usize::MAX, workspace);
self.total += u64::from(size).pow(2);
let id = self.free_ids.pop().unwrap_or_else(|| {
self.members.push(Vec::new());
self.members.len() - 1
});
debug_assert!(self.members[id].is_empty());
for &(v, separated, constant) in &workspace.visited {
let vi = CoordIndex(v as usize);
let gain = 2 * size * u32::from(separated) - constant;
self.component[vi] = id as u16;
self.gain[vi] = gain;
gains[vi] += u64::from(gain);
self.members[id].push(v);
}
}
fn set_cut(
&mut self,
v: usize,
cut: bool,
gains: &mut Map2d<u64>,
workspace: &mut Workspace,
) {
let vi = CoordIndex(v);
if !self.dry[vi] {
return;
}
debug_assert_eq!(self.order[vi] == u16::MAX, !cut);
let mut roots = [0; 4];
let mut count = 0;
if cut {
roots[0] = self.component[vi] as usize;
count = 1;
} else {
for next in [v - N, v + 1, v + N, v - 1] {
let ni = CoordIndex(next);
if self.order[ni] != u16::MAX {
let root = self.component[ni] as usize;
if !roots[..count].contains(&root) {
roots[count] = root;
count += 1;
}
}
}
}
workspace.affected.clear();
for &root in &roots[..count] {
self.total -= (self.members[root].len() as u64).pow(2);
for &u in &self.members[root] {
let ui = CoordIndex(u as usize);
gains[ui] -= u64::from(self.gain[ui]);
self.gain[ui] = 0;
self.order[ui] = 0;
workspace.affected.push(u);
}
self.members[root].clear();
self.free_ids.push(root);
}
self.order[vi] = if cut { u16::MAX } else { 0 };
if !cut {
workspace.affected.push(v as u16);
}
for i in 0..workspace.affected.len() {
let root = workspace.affected[i] as usize;
if self.order[CoordIndex(root)] == 0 {
self.rebuild_component(root, gains, workspace);
}
}
}
fn restoration_cost(&self, v: usize) -> u64 {
if !self.dry[CoordIndex(v)] {
return 0;
}
let mut roots = [0; 4];
let mut count = 0;
let mut size = 1_u64;
let mut squares = 0;
for next in [v - N, v + 1, v + N, v - 1] {
let ni = CoordIndex(next);
if self.order[ni] == u16::MAX {
continue;
}
let root = self.component[ni] as usize;
if !roots[..count].contains(&root) {
roots[count] = root;
count += 1;
let len = self.members[root].len() as u64;
size += len;
squares += len * len;
}
}
size * size - squares
}
fn exchange_delta(
&mut self,
removed: &[Coord],
added: &[Coord],
workspace: &mut Workspace,
) -> i64 {
let mut roots = [usize::MAX; 16];
let mut count = 0;
for &p in removed {
let v = index(p).0;
if !self.dry[CoordIndex(v)] {
continue;
}
for next in [v - N, v + 1, v + N, v - 1] {
let ni = CoordIndex(next);
if self.order[ni] == u16::MAX {
continue;
}
let root = self.component[ni] as usize;
if !roots[..count].contains(&root) {
roots[count] = root;
count += 1;
}
}
}
for &p in added {
let vi = index(p);
if !self.dry[vi] {
continue;
}
let root = self.component[vi] as usize;
if !roots[..count].contains(&root) {
roots[count] = root;
count += 1;
}
}
workspace.affected.clear();
let mut old = 0_i64;
for &root in &roots[..count] {
old += (self.members[root].len() as i64).pow(2);
for &v in &self.members[root] {
workspace.affected.push(v);
self.order[CoordIndex(v as usize)] = 0;
}
}
for &p in removed {
let vi = index(p);
if self.dry[vi] {
self.order[vi] = 0;
workspace.affected.push(vi.0 as u16);
}
}
for &p in added {
self.order[index(p)] = u16::MAX;
}
let mut new = 0_i64;
for i in 0..workspace.affected.len() {
let root = workspace.affected[i] as usize;
if self.order[CoordIndex(root)] != 0 {
continue;
}
workspace.queue.clear();
workspace.queue.push(root as u16);
self.order[CoordIndex(root)] = 1;
let mut head = 0;
while head < workspace.queue.len() {
let v = workspace.queue[head] as usize;
head += 1;
for next in [v - N, v + 1, v + N, v - 1] {
let ni = CoordIndex(next);
if self.order[ni] == 0 {
self.order[ni] = 1;
workspace.queue.push(next as u16);
}
}
}
new += (workspace.queue.len() as i64).pow(2);
}
// 未採用の候補なので、乾燥頂点の有無だけ元に戻す。DFS順序の値は再計算時に初期化する。
for &p in removed {
self.order[index(p)] = u16::MAX;
}
for &p in added {
if self.dry[index(p)] {
self.order[index(p)] = 1;
}
}
new - old
}
}
pub(super) struct IncrementalEvaluator {
samples: Vec<Sample>,
workspace: Workspace,
pub(super) gains: Map2d<u64>,
pub(super) total: u64,
}
impl IncrementalEvaluator {
pub(super) fn new(dry: &[Map2d<bool>], cuts: &[Coord]) -> Self {
let mut gains = Map2d::with_default(N);
let mut workspace = Workspace::new();
let samples: Vec<_> = dry
.iter()
.map(|d| Sample::new(d, cuts, &mut gains, &mut workspace))
.collect();
let total = samples.iter().map(|s| s.total).sum();
Self {
samples,
workspace,
gains,
total,
}
}
pub(super) fn extend(&mut self, dry: &[Map2d<bool>], cuts: &[Coord]) {
for d in dry {
let sample = Sample::new(d, cuts, &mut self.gains, &mut self.workspace);
self.total += sample.total;
self.samples.push(sample);
}
}
pub(super) fn set_cut(&mut self, p: Coord, cut: bool) {
self.total = 0;
let v = index(p).0;
for sample in &mut self.samples {
sample.set_cut(v, cut, &mut self.gains, &mut self.workspace);
self.total += sample.total;
}
}
/// 伐採済みの頂点を1本戻したときの二乗和の増加。
#[cfg(test)]
pub(super) fn restoration_cost(&self, p: Coord) -> u64 {
let v = index(p).0;
self.samples
.iter()
.map(|sample| sample.restoration_cost(v))
.sum()
}
/// サンプルごとのデータをまとめて参照し、候補全体の復元コストを求める。
pub(super) fn restoration_costs(&self, cuts: &[Coord]) -> Vec<u64> {
let indices: Vec<_> = cuts.iter().map(|&p| index(p).0).collect();
let mut costs = vec![0; cuts.len()];
for sample in &self.samples {
for (&v, cost) in indices.iter().zip(&mut costs) {
*cost += sample.restoration_cost(v);
}
}
costs
}
/// 状態を採用せずに、2本以内の交換後の二乗和を求める。
pub(super) fn exchange_total(&mut self, removed: &[Coord], added: &[Coord]) -> u64 {
assert!(removed.len() <= 2 && added.len() <= 2);
let delta: i64 = self
.samples
.iter_mut()
.map(|sample| sample.exchange_delta(removed, added, &mut self.workspace))
.sum();
(self.total as i64 + delta) as u64
}
}
}
mod refinement {
//! 固定サンプル上の頂点削除利得をまとめて求め、伐採1本の交換で仕上げる。
use super::evaluation::Evaluator;
use super::incremental::IncrementalEvaluator;
use crate::library::{
diagnostics::Timer,
grid::{Coord, CoordIndex, Map2d},
random::Random,
};
use crate::problem::{Input, Output};
use std::time::Duration;
const N: usize = Input::N + 2;
const RESORT_INTERVAL: usize = 50;
const INITIAL_SAMPLES: usize = 128;
const FINAL_SAMPLE_LEAD_MS: u64 = 170;
fn index(p: Coord) -> CoordIndex {
CoordIndex((p.row() + 1) * N + p.col() + 1)
}
#[cfg(test)]
struct Workspace {
order: Map2d<u16>,
// (頂点, 1 + 切り離す子部分木の総サイズ, その二乗 + 子部分木サイズの二乗和)
visited: Vec<(u16, u16, u32)>,
clock: u16,
}
#[cfg(test)]
impl Workspace {
fn new() -> Self {
Self {
order: Map2d::from_fn(|_| u16::MAX, N),
visited: Vec::with_capacity(Input::N * Input::N),
clock: 0,
}
}
fn dfs(&mut self, v: usize, parent: usize) -> (u16, u32) {
self.clock += 1;
let order = self.clock;
self.order[CoordIndex(v)] = order;
let mut low = order;
let mut size = 1;
let mut separated_size = 1;
let mut separated_square = 0;
for next in [v - N, v + 1, v + N, v - 1] {
let next_order = self.order[CoordIndex(next)];
if next_order == 0 {
let (child_low, child_size) = self.dfs(next, v);
size += child_size;
low = low.min(child_low);
if child_low >= order {
separated_size += child_size;
separated_square += child_size * child_size;
}
} else if next != parent {
// 非乾燥・伐採マスと外周はMAXなので、lowを更新しない。
low = low.min(next_order);
}
}
self.visited.push((
v as u16,
separated_size as u16,
separated_size * separated_size + separated_square,
));
(low, size)
}
/// 現在の二乗和と、各頂点を追加で伐採したときの二乗和の減少量。
fn analyze(
&mut self,
samples: &[Map2d<bool>],
cuts: &[Coord],
gains: &mut Map2d<u64>,
) -> u64 {
for row in 0..N {
gains[row].fill(0);
}
let mut total = 0;
for sample in samples {
for row in 0..Input::N {
for col in 0..Input::N {
self.order[row + 1][col + 1] =
if sample[row][col] { 0 } else { u16::MAX };
}
}
for &p in cuts {
self.order[index(p)] = u16::MAX;
}
self.clock = 0;
self.visited.clear();
for start in N..N * (Input::N + 1) {
let si = CoordIndex(start);
if self.order[si] != 0 {
continue;
}
let begin = self.visited.len();
let (_, component_size) = self.dfs(start, usize::MAX);
total += u64::from(component_size).pow(2);
for &(v, separated, constant) in &self.visited[begin..] {
// S² - (S - separated)² - 子部分木の二乗和を展開した形。
gains[CoordIndex(v as usize)] +=
2 * u64::from(component_size) * u64::from(separated)
- u64::from(constant);
}
}
}
total
}
}
pub(super) fn refine(
input: &Input,
mut output: Output,
timer: &Timer,
deadline: Duration,
rng: &mut Random,
evaluator: &Evaluator,
) -> Output {
if timer.elapsed() >= deadline {
return output;
}
let vertices: Vec<_> = (0..Input::N)
.flat_map(|row| (0..Input::N).map(move |col| Coord::new(row, col)))
.filter(|&p| input.d[p] > 0)
.collect();
if output.cuts.len() == vertices.len() {
return output;
}
let max_samples = evaluator.samples().len();
let mut sample_count = INITIAL_SAMPLES.min(max_samples);
let mut objective =
IncrementalEvaluator::new(&evaluator.samples()[..sample_count], &output.cuts);
let mut current = objective.total;
let mut cut = Map2d::with_default(Input::N);
for &p in &output.cuts {
cut[p] = true;
}
let mut steps = 0;
let mut improvements = 0;
let mut stalled = 0;
let mut permutation: Vec<_> = (0..output.cuts.len()).collect();
let mut costs: Vec<_> = objective.restoration_costs(&output.cuts);
permutation.sort_by_key(|&i| costs[i]);
let mut cursor = 0;
while timer.elapsed() < deadline {
let elapsed = timer.elapsed();
let target = if elapsed
>= deadline.saturating_sub(Duration::from_millis(FINAL_SAMPLE_LEAD_MS))
{
max_samples
} else if elapsed
>= deadline
.saturating_sub(Duration::from_millis(FINAL_SAMPLE_LEAD_MS + 200))
{
(INITIAL_SAMPLES * 2).min(max_samples)
} else if stalled >= Input::K {
(sample_count * 2).min(max_samples)
} else {
sample_count
};
if target > sample_count {
objective.extend(&evaluator.samples()[sample_count..target], &output.cuts);
sample_count = target;
current = objective.total;
stalled = 0;
cursor = 0;
costs = objective.restoration_costs(&output.cuts);
permutation.sort_by_key(|&i| costs[i]);
} else if stalled >= Input::K && sample_count == max_samples {
break;
}
let removed = if output.cuts.len() == Input::K {
let slot = permutation[cursor % permutation.len()];
let p = output.cuts.swap_remove(slot);
cut[p] = false;
objective.set_cut(p, false);
Some((slot, p))
} else {
None
};
let total = objective.total;
let mut chosen = removed.map(|(_, p)| p).unwrap_or(vertices[0]);
for &p in &vertices {
if !cut[p]
&& (cut[chosen]
|| objective.gains[index(p)] > objective.gains[index(chosen)])
{
chosen = p;
}
}
let next = total - objective.gains[index(chosen)];
debug_assert!(next <= current);
if next < current {
improvements += 1;
stalled = 0;
} else {
stalled += 1;
}
current = next;
cut[chosen] = true;
objective.set_cut(chosen, true);
debug_assert_eq!(objective.total, next);
output.cuts.push(chosen);
if output.cuts.len() == vertices.len() {
break;
}
if let Some((slot, _)) = removed {
let last = output.cuts.len() - 1;
output.cuts.swap(slot, last);
} else {
permutation.push(permutation.len());
}
steps += 1;
cursor += 1;
if !permutation.is_empty()
&& (cursor % permutation.len() == 0 || cursor % RESORT_INTERVAL == 0)
{
costs = objective.restoration_costs(&output.cuts);
let first = cursor % permutation.len();
permutation[first..].sort_by_key(|&i| costs[i]);
}
}
let mut flips = 0;
let mut flip_trials = 0;
while timer.elapsed() < deadline {
let mut blocks = Vec::new();
for row in 0..Input::N - 1 {
for col in 0..Input::N - 1 {
let block = [
Coord::new(row, col),
Coord::new(row + 1, col),
Coord::new(row, col + 1),
Coord::new(row + 1, col + 1),
];
if block.iter().filter(|&&p| cut[p]).count() == 2
&& block.iter().all(|&p| input.d[p] > 0)
{
blocks.push(block);
}
}
}
for i in (1..blocks.len()).rev() {
let j = rng.index(i + 1);
blocks.swap(i, j);
}
let mut improved = false;
for block in blocks {
if timer.elapsed() >= deadline {
break;
}
let removed: Vec<_> = block.iter().copied().filter(|&p| cut[p]).collect();
if removed.len() != 2 {
continue;
}
let added: Vec<_> = block.iter().copied().filter(|&p| !cut[p]).collect();
let next = objective.exchange_total(&removed, &added);
flip_trials += 1;
if next < current {
for &p in &removed {
cut[p] = false;
objective.set_cut(p, false);
}
for &p in &added {
cut[p] = true;
objective.set_cut(p, true);
}
debug_assert_eq!(next, objective.total);
current = next;
flips += 1;
improved = true;
}
}
if !improved {
break;
}
}
output.cuts = vertices.iter().copied().filter(|&p| cut[p]).collect();
eprintln!("local_flip_trials={flip_trials} flips={flips} samples={sample_count}");
eprintln!(
"refinement_steps={steps} improvements={improvements} refinement_square_sum={current}"
);
output
}
#[cfg(test)]
mod tests {
use super::*;
use crate::library::grid::ADJACENTS;
#[test]
fn empty_or_sparse_forest_does_not_duplicate_cuts() {
for trees in [0, 2] {
let input = Input {
d: Map2d::from_fn(
|p| {
if p.row() == 0 && p.col() < trees {
9
} else {
0
}
},
Input::N,
),
};
let evaluator = Evaluator::with_samples(&input, 512);
let output = refine(
&input,
Output::default(),
&Timer::default(),
Duration::from_secs(1),
&mut Random::new(42),
&evaluator,
);
let expected: Vec<_> = (0..trees).map(|col| Coord::new(0, col)).collect();
assert_eq!(output.cuts, expected);
}
}
#[test]
fn incremental_updates_match_full_recalculation() {
let mut rng = Random::new(12345);
let samples: Vec<_> = [2, 7, 10]
.into_iter()
.map(|density| {
Map2d::from_fn(
|p| p.row() < 10 && p.col() < 10 && rng.index(10) < density,
Input::N,
)
})
.collect();
let mut cuts = vec![Coord::new(0, 0), Coord::new(5, 5)];
let mut objective = IncrementalEvaluator::new(&samples[..2], &cuts);
let mut workspace = Workspace::new();
let mut gains = Map2d::with_default(N);
for iteration in 0..300 {
if iteration == 100 {
objective.extend(&samples[2..], &cuts);
}
let active_samples = &samples[..if iteration < 100 { 2 } else { 3 }];
let total = workspace.analyze(active_samples, &cuts, &mut gains);
assert_eq!(objective.total, total);
assert!(objective.gains.iter().eq(gains.iter()));
if cuts.len() >= 2 {
let removed = &cuts[..2];
let added: Vec<_> = (0..121)
.map(|i| Coord::new(i / 11, i % 11))
.filter(|p| !cuts.contains(p))
.take(2)
.collect();
if added.len() == 2 {
let mut candidate = cuts[2..].to_vec();
candidate.extend_from_slice(&added);
let expected: u64 = active_samples
.iter()
.map(|sample| brute_sum(sample, &candidate))
.sum();
assert_eq!(objective.exchange_total(removed, &added), expected);
assert_eq!(objective.total, total);
assert!(objective.gains.iter().eq(gains.iter()));
}
}
assert_eq!(
objective.restoration_costs(&cuts),
cuts.iter()
.map(|&p| objective.restoration_cost(p))
.collect::<Vec<_>>()
);
let p = Coord::new(rng.index(11), rng.index(11));
if let Some(i) = cuts.iter().position(|&q| q == p) {
cuts.swap_remove(i);
let cost = objective.restoration_cost(p);
let old = objective.total;
objective.set_cut(p, false);
assert_eq!(objective.total, old + cost);
} else {
cuts.push(p);
objective.set_cut(p, true);
}
}
}
fn brute_sum(sample: &Map2d<bool>, cuts: &[Coord]) -> u64 {
let mut seen = Map2d::from_fn(|p| !sample[p], Input::N);
for &p in cuts {
seen[p] = true;
}
let mut total = 0;
for row in 0..Input::N {
for col in 0..Input::N {
let start = Coord::new(row, col);
if seen[start] {
continue;
}
seen[start] = true;
let mut queue = vec![start];
let mut head = 0;
while head < queue.len() {
let p = queue[head];
head += 1;
for diff in ADJACENTS {
let q = p + diff;
if q.in_map(Input::N) && !seen[q] {
seen[q] = true;
queue.push(q);
}
}
}
total += (queue.len() as u64).pow(2);
}
}
total
}
#[test]
fn all_vertex_gains_match_brute_force() {
let mut rng = Random::new(2026);
let mut workspace = Workspace::new();
let mut gains = Map2d::with_default(N);
for density in [3, 7, 10] {
let sample = Map2d::from_fn(
|p| p.row() < 10 && p.col() < 10 && rng.index(10) < density,
Input::N,
);
let mut cuts = vec![Coord::new(5, 5), Coord::new(0, 0)];
let total =
workspace.analyze(std::slice::from_ref(&sample), &cuts, &mut gains);
assert_eq!(total, brute_sum(&sample, &cuts));
for row in 0..11 {
for col in 0..11 {
let p = Coord::new(row, col);
cuts.push(p);
assert_eq!(
gains[index(p)],
total - brute_sum(&sample, &cuts),
"{p:?}, density={density}"
);
cuts.pop();
}
}
}
}
}
}
use std::{collections::HashMap, time::Duration};
use crate::library::{
bitboard::BitBoard,
diagnostics::Timer,
grid::{ADJACENTS, Coord, CoordIndex, Map2d},
maxflow::MaxFlow,
random::Random,
};
use crate::problem::{Input, Output};
use evaluation::Evaluator;
const D_THRESHOLD: u8 = 5;
// 乾燥度を優先しつつ、少ない本数で分断できるカットを選ぶ。
const VERTEX_BASE_COST: i64 = 20;
const PAIR_TRIALS: usize = 1;
const RADIUS: usize = 12;
const MIN_RADIUS: usize = RADIUS;
const MAX_RADIUS: usize = RADIUS;
#[cfg(test)]
const TERMINAL_RADIUS: usize = 12;
const REMOVE_LOWER: usize = 1;
const REMOVE_UPPER: usize = 2;
// 未使用の伐採枠は仕上げの1本交換で埋められるため、構築途中でもその余地を評価する。
const UNUSED_CUT_BONUS: i64 = 15000;
const TIME_LIMIT: Duration = Duration::from_millis(1980);
const REFINEMENT_SAMPLES: usize = 1024;
const CONSTRUCTION_LIMIT: Duration = Duration::from_millis(1200);
fn components(input: &Input, cut: &Map2d<bool>) -> Vec<Vec<Coord>> {
// 外周を非対象マスにして境界判定を省く。成分・頂点の順序は元のBFSと揃える。
const PAD: usize = Input::N + 2;
let mut active = Map2d::with_default(PAD);
for row in 0..Input::N {
for col in 0..Input::N {
active[row + 1][col + 1] = input.d[row][col] >= D_THRESHOLD && !cut[row][col];
}
}
let mut result = Vec::new();
let mut queue = Vec::with_capacity(Input::N * Input::N);
for start in PAD..PAD * (Input::N + 1) {
if !active[CoordIndex(start)] {
continue;
}
active[CoordIndex(start)] = false;
queue.clear();
queue.push(start);
let mut head = 0;
while head < queue.len() {
let v = queue[head];
head += 1;
for next in [v - PAD, v + 1, v + PAD, v - 1] {
if active[CoordIndex(next)] {
active[CoordIndex(next)] = false;
queue.push(next);
}
}
}
result.push(
queue
.iter()
.map(|&v| Coord::new(v / PAD - 1, v % PAD - 1))
.collect(),
);
}
result
}
/// 乾燥度を重視した成分の重みを累積し、抽選ごとの合計計算を省く。
fn component_weights(input: &Input, components: &[Vec<Coord>]) -> Vec<usize> {
components
.iter()
.map(|component| {
let mass: usize = component
.iter()
.map(|&p| usize::from(input.d[p] - 4).pow(2))
.sum();
mass * mass
})
.scan(0, |sum, weight| {
*sum += weight;
Some(*sum)
})
.collect()
}
fn choose_component<'a>(
components: &'a [Vec<Coord>],
cumulative_weights: &[usize],
rng: &mut Random,
) -> Option<&'a [Coord]> {
let total = cumulative_weights.last().copied().unwrap_or(0);
if total == 0 {
return None;
}
let draw = rng.index(total);
let i = cumulative_weights.partition_point(|&prefix| prefix <= draw);
Some(&components[i])
}
fn choose_pair(vertices: &[Coord], rng: &mut Random) -> (Coord, Coord) {
let mut best = (vertices[0], vertices[0]);
for _ in 0..PAIR_TRIALS {
let pair = (
vertices[rng.index(vertices.len())],
vertices[rng.index(vertices.len())],
);
if pair.0.dist(&pair.1) > best.0.dist(&best.1) {
best = pair;
}
}
best
}
fn minimum_vertex_cut(
input: &Input,
vertices: &[Coord],
p1: Coord,
p2: Coord,
radius: usize,
) -> Vec<Coord> {
let len = vertices.len();
let mut indices = Map2d::from_fn(|_| usize::MAX, Input::N);
let mut board = BitBoard::<{ Input::N }>::default();
for (i, &p) in vertices.iter().enumerate() {
indices[p] = i;
board.insert(p);
}
let source_vertices = board.within_distance(p1, radius);
let sink_vertices = board.within_distance(p2, radius);
let source = 2 * len;
let sink = source + 1;
let mut graph = MaxFlow::new(sink + 1);
// 全頂点の有限コストの総和より大きく、隣接辺はカットされない。
let infinity = (len as i64) * VERTEX_BASE_COST + 1;
for (i, &p) in vertices.iter().enumerate() {
let vin = 2 * i;
let vout = vin + 1;
graph.add_edge(vin, vout, VERTEX_BASE_COST - i64::from(input.d[p]));
// source/sink内の頂点も有限コストのvin→voutを通るので伐採可能。
if source_vertices.contains(p) {
graph.add_edge(source, vin, infinity);
}
if sink_vertices.contains(p) {
graph.add_edge(vout, sink, infinity);
}
for diff in ADJACENTS {
let next = p + diff;
if next.in_map(Input::N) && indices[next] != usize::MAX {
graph.add_edge(vout, 2 * indices[next], infinity);
}
}
}
graph.flow(source, sink);
let reachable = graph.min_cut(source);
vertices
.iter()
.enumerate()
.filter_map(|(i, &p)| (reachable[2 * i] && !reachable[2 * i + 1]).then_some(p))
.collect()
}
#[derive(Clone)]
struct State {
cut: Map2d<bool>,
batches: Vec<Vec<Coord>>,
count: usize,
}
impl State {
fn new() -> Self {
Self {
cut: Map2d::with_default(Input::N),
batches: Vec::new(),
count: 0,
}
}
fn add(&mut self, batch: Vec<Coord>) {
for &p in &batch {
debug_assert!(!self.cut[p]);
self.cut[p] = true;
}
self.count += batch.len();
self.batches.push(batch);
}
fn remove_random_batches(&mut self, rng: &mut Random) {
let number = REMOVE_LOWER + rng.index(REMOVE_UPPER - REMOVE_LOWER + 1);
// 残っている履歴から毎回一様に選び、重複なしでカットを取り消す。
for _ in 0..number.min(self.batches.len()) {
let index = rng.index(self.batches.len());
let removed = self.batches.swap_remove(index);
self.count -= removed.len();
for p in removed {
self.cut[p] = false;
}
}
}
fn output(&self) -> Output {
Output {
cuts: self.batches.iter().flatten().copied().collect(),
}
}
fn key(&self) -> [u64; Input::N] {
let mut result = [0; Input::N];
for &p in self.batches.iter().flatten() {
result[p.row()] |= 1_u64 << p.col();
}
result
}
}
fn complete(input: &Input, state: &mut State, rng: &mut Random, timer: &Timer) {
let mut current = components(input, &state.cut);
let mut weights = component_weights(input, ¤t);
// 上限超過・再抽選の上限・構築フェーズの締切で終了する。
let mut failures = 0;
loop {
failures += 1;
if failures > 1024 || timer.elapsed() >= CONSTRUCTION_LIMIT {
break;
}
let Some(vertices) = choose_component(¤t, &weights, rng) else {
continue;
};
let radius = MIN_RADIUS + rng.index(MAX_RADIUS - MIN_RADIUS + 1);
let (p1, p2) = choose_pair(vertices, rng);
if p1.dist(&p2) < 2 * radius {
continue;
}
let added = minimum_vertex_cut(input, vertices, p1, p2, radius);
if state.count + added.len() > Input::K {
break;
}
state.add(added);
failures = 0;
current = components(input, &state.cut);
weights = component_weights(input, ¤t);
}
}
fn complete_greedy(
input: &Input,
state: &mut State,
rng: &mut Random,
evaluator: &Evaluator,
timer: &Timer,
) {
let mut value = evaluator.total(&state.cut);
while state.count < Input::K {
let current = components(input, &state.cut);
let weights = component_weights(input, ¤t);
let mut best = None;
let mut best_efficiency = 0.0;
let mut candidates = 0;
for _ in 0..1024 {
if timer.elapsed() >= CONSTRUCTION_LIMIT {
break;
}
let Some(vertices) = choose_component(¤t, &weights, rng) else {
break;
};
let (p1, p2) = choose_pair(vertices, rng);
let radius = RADIUS;
if p1.dist(&p2) < 2 * radius {
continue;
}
let batch = minimum_vertex_cut(input, vertices, p1, p2, radius);
if batch.is_empty() || state.count + batch.len() > Input::K {
continue;
}
for &p in &batch {
state.cut[p] = true;
}
let next = evaluator.total(&state.cut);
for &p in &batch {
state.cut[p] = false;
}
let efficiency = (value - next) as f64 / batch.len() as f64;
if efficiency > best_efficiency {
best_efficiency = efficiency;
best = Some((batch, next));
}
candidates += 1;
if candidates == 8 {
break;
}
}
let Some((batch, next)) = best else {
break;
};
state.add(batch);
value = next;
}
}
fn accept_candidate(
state: &mut State,
score: &mut i64,
candidate: State,
candidate_score: i64,
allow_worse: bool,
) {
// 却下時は現在解を残し、受理時だけ盤面・本数・カット履歴をまとめて置き換える。
if candidate_score >= *score || allow_worse {
*state = candidate;
*score = candidate_score;
}
}
struct Archived {
state: State,
score: i64,
}
fn remember(archive: &mut Vec<Archived>, candidate: &State, score: i64) {
if archive.len() == 16 && archive.iter().all(|e| e.score >= score) {
return;
}
for entry in archive.iter_mut() {
let common = candidate
.batches
.iter()
.flatten()
.filter(|&&p| entry.state.cut[p])
.count();
if candidate.count + entry.state.count - 2 * common < 12 {
if score > entry.score {
*entry = Archived {
state: candidate.clone(),
score,
};
}
return;
}
}
if archive.len() < 16 {
archive.push(Archived {
state: candidate.clone(),
score,
});
} else {
let i = archive
.iter()
.enumerate()
.min_by_key(|(_, e)| e.score)
.unwrap()
.0;
if score > archive[i].score {
archive[i] = Archived {
state: candidate.clone(),
score,
};
}
}
}
pub fn solve(input: &Input, timer: &Timer) -> Output {
let evaluator = Evaluator::new(input);
let mut rng = Random::new(42);
let mut best = Output::default();
let mut archive = Vec::new();
let mut best_score = i64::MIN;
let mut trials = 0;
let mut state = State::new();
let mut current_score = i64::MIN;
let mut score_cache = HashMap::new();
let mut cache_hits = 0;
loop {
let mut candidate = state.clone();
if trials > 0 {
candidate.remove_random_batches(&mut rng);
}
if trials == 0 {
complete_greedy(input, &mut candidate, &mut rng, &evaluator, timer);
eprintln!(
"initial_elapsed_ms={:.3}",
timer.elapsed().as_secs_f64() * 1000.0
);
} else {
complete(input, &mut candidate, &mut rng, timer);
}
let key = candidate.key();
let raw_score = if let Some(&cached) = score_cache.get(&key) {
cache_hits += 1;
cached
} else {
let score = evaluator.score(&candidate.cut);
score_cache.insert(key, score);
score
};
let score = raw_score + UNUSED_CUT_BONUS * (Input::K - candidate.count) as i64;
trials += 1;
remember(&mut archive, &candidate, score);
if score > best_score {
best_score = score;
best = candidate.output();
}
let progress =
(timer.elapsed().as_secs_f64() / CONSTRUCTION_LIMIT.as_secs_f64()).min(1.0);
let temperature = 10000_f64 * 0.1_f64.powf(progress);
let allow_worse = score < current_score
&& ((rng.next_u64() >> 11) as f64 / (1_u64 << 53) as f64)
< ((score - current_score) as f64 / temperature).exp();
accept_candidate(
&mut state,
&mut current_score,
candidate,
score,
allow_worse,
);
// 評価完了後だけ時間を判定する。探索用乱数の状態は巻き戻さない。
if timer.elapsed() > CONSTRUCTION_LIMIT {
break;
}
}
let validation = Evaluator::with_samples(input, REFINEMENT_SAMPLES);
if let Some(entry) = archive.iter().max_by_key(|entry| {
validation.score(&entry.state.cut)
+ UNUSED_CUT_BONUS * (Input::K - entry.state.count) as i64
}) {
best = entry.state.output();
}
eprintln!("construction_cuts={}", best.cuts.len());
eprintln!(
"cache_hits={cache_hits} cache_entries={}",
score_cache.len()
);
best = refinement::refine(input, best, timer, TIME_LIMIT, &mut rng, &validation);
eprintln!(
"trials={trials} construction_objective={best_score} cuts={} elapsed_ms={:.3}",
best.cuts.len(),
timer.elapsed().as_secs_f64() * 1000.0
);
best
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn padded_components_preserve_component_and_vertex_order() {
let mut rng = Random::new(987);
for density in [0, 5, 8, 10] {
for _ in 0..8 {
let input = Input {
d: Map2d::from_fn(
|_| {
if rng.index(10) < density {
5 + rng.index(5) as u8
} else {
0
}
},
Input::N,
),
};
let cut = Map2d::from_fn(|_| rng.index(7) == 0, Input::N);
assert_eq!(components(&input, &cut), components_reference(&input, &cut));
}
}
}
#[test]
fn component_weights_accumulate_squared_dryness_mass() {
let mut input = Input {
d: Map2d::with_default(Input::N),
};
input.d[Coord::new(0, 0)] = 5;
input.d[Coord::new(0, 1)] = 6;
input.d[Coord::new(3, 3)] = 9;
input.d[Coord::new(6, 6)] = 7;
input.d[Coord::new(0, 2)] = 4;
let cut = Map2d::with_default(Input::N);
let components = components(&input, &cut);
// 成分の重みは(1+4)^2、25^2、9^2。その累積和を返す。
assert_eq!(component_weights(&input, &components), [25, 650, 731]);
}
#[test]
fn cumulative_selection_matches_weight_subtraction() {
let components: Vec<_> = (0..7).map(|row| vec![Coord::new(row, 0)]).collect();
let raw = [0, 11, 3, 0, 19, 2, 0];
let cumulative: Vec<_> = raw
.iter()
.scan(0, |sum, weight| {
*sum += weight;
Some(*sum)
})
.collect();
let mut actual_rng = Random::new(42);
let mut expected_rng = Random::new(42);
for _ in 0..10000 {
let mut draw = expected_rng.index(raw.iter().sum());
let expected = raw
.iter()
.position(|&weight| {
if draw < weight {
true
} else {
draw -= weight;
false
}
})
.unwrap();
assert_eq!(
choose_component(&components, &cumulative, &mut actual_rng),
Some(components[expected].as_slice())
);
}
assert!(choose_component(&[], &[], &mut actual_rng).is_none());
}
#[test]
fn component_selection_obeys_zero_weights() {
let components = vec![vec![Coord::new(0, 0)], vec![Coord::new(1, 1)]];
let mut rng = Random::new(42);
for _ in 0..100 {
assert_eq!(
choose_component(&components, &[0, 1], &mut rng),
Some(components[1].as_slice())
);
}
assert!(choose_component(&components, &[0, 0], &mut rng).is_none());
}
#[test]
fn score_cache_key_ignores_batch_order() {
let mut a = State::new();
a.add(vec![Coord::new(0, 0), Coord::new(49, 49)]);
a.add(vec![Coord::new(12, 24)]);
let mut b = State::new();
b.add(vec![Coord::new(12, 24), Coord::new(49, 49)]);
b.add(vec![Coord::new(0, 0)]);
assert_eq!(a.key(), b.key());
b.add(vec![Coord::new(1, 1)]);
assert_ne!(a.key(), b.key());
}
#[test]
fn archive_keeps_high_scores_and_replaces_similar_candidate() {
let mut archive = Vec::new();
for row in 0..20 {
let mut state = State::new();
state.add((0..7).map(|col| Coord::new(row, col)).collect());
remember(&mut archive, &state, row as i64);
}
assert_eq!(archive.len(), 16);
assert!(archive.iter().all(|entry| entry.score >= 4));
let mut candidate = State::new();
candidate.add((0..7).map(|col| Coord::new(19, col)).collect());
remember(&mut archive, &candidate, 100);
assert_eq!(archive.len(), 16);
assert_eq!(archive.iter().filter(|entry| entry.score == 100).count(), 1);
}
#[test]
fn annealing_can_accept_a_worsening_candidate() {
let mut state = State::new();
state.add(vec![Coord::new(0, 0)]);
let mut score = 100;
let mut candidate = State::new();
candidate.add(vec![Coord::new(1, 1)]);
accept_candidate(&mut state, &mut score, candidate, 99, true);
assert_eq!(score, 99);
assert_eq!(state.output().cuts, vec![Coord::new(1, 1)]);
}
fn line_input_with_existing_cuts(count: usize) -> (Input, State) {
let mut input = Input {
d: Map2d::with_default(Input::N),
};
for col in 0..Input::N {
input.d[Coord::new(0, col)] = 5;
}
let existing: Vec<_> = (0..count)
.map(|i| Coord::new(2 + i / Input::N, i % Input::N))
.collect();
for &p in &existing {
input.d[p] = 1;
}
let mut state = State::new();
state.add(existing);
(input, state)
}
#[test]
fn continues_after_95_and_adds_a_cut_that_reaches_k() {
let (input, mut state) = line_input_with_existing_cuts(Input::K - 1);
complete(&input, &mut state, &mut Random::new(42), &Timer::default());
assert_eq!(state.count, Input::K);
assert_eq!(state.batches.len(), 2);
assert_eq!(state.batches[1].len(), 1);
assert_eq!(state.cut[0].iter().filter(|&&cut| cut).count(), 1);
}
#[test]
fn overflowing_cut_ends_trial_without_changing_state() {
let (input, mut state) = line_input_with_existing_cuts(Input::K);
let saved = state.clone();
complete(&input, &mut state, &mut Random::new(42), &Timer::default());
assert_eq!(state.count, saved.count);
assert_eq!(state.batches, saved.batches);
assert!(state.cut.iter().eq(saved.cut.iter()));
}
#[test]
fn worsening_candidate_preserves_state_and_history() {
let mut state = State::new();
state.add(vec![Coord::new(0, 0), Coord::new(0, 1)]);
state.add(vec![Coord::new(1, 0)]);
let saved = state.clone();
let mut score = 100;
let mut candidate = State::new();
candidate.add(vec![Coord::new(2, 2)]);
accept_candidate(&mut state, &mut score, candidate, 99, false);
assert_eq!(score, 100);
assert_eq!(state.batches, saved.batches);
assert_eq!(state.count, saved.count);
assert!(state.cut.iter().eq(saved.cut.iter()));
}
#[test]
fn equal_and_better_candidates_are_accepted() {
let mut state = State::new();
state.add(vec![Coord::new(0, 0)]);
let mut score = 100;
for (col, candidate_score) in [(1, 100), (2, 101)] {
let mut candidate = State::new();
candidate.add(vec![Coord::new(0, col)]);
accept_candidate(&mut state, &mut score, candidate, candidate_score, false);
assert_eq!(score, candidate_score);
assert_eq!(state.output().cuts, [Coord::new(0, col)]);
assert_eq!(state.count, 1);
assert_eq!(state.cut.iter().filter(|&&value| value).count(), 1);
}
}
#[test]
fn components_exclude_low_dryness_and_cut_cells() {
let mut input = Input {
d: Map2d::with_default(Input::N),
};
for (col, d) in [5, 9, 4, 5, 5, 9].into_iter().enumerate() {
input.d[Coord::new(0, col)] = d;
}
let mut cut = Map2d::with_default(Input::N);
cut[Coord::new(0, 4)] = true;
let actual = components(&input, &cut);
assert_eq!(actual.iter().map(Vec::len).collect::<Vec<_>>(), [2, 1, 1]);
assert_eq!(actual[0], [Coord::new(0, 0), Coord::new(0, 1)]);
}
#[test]
fn terminal_sets_follow_component_detour() {
let mut input = Input {
d: Map2d::from_fn(|_| 5, Input::N),
};
let mut vertices = Vec::new();
for row in 0..16 {
vertices.push(Coord::new(row, 0));
vertices.push(Coord::new(row, 2));
}
let bridge = Coord::new(15, 1);
vertices.push(bridge);
input.d[bridge] = 9;
// 2点のマンハッタン距離は2だが、成分内のBFS距離は32。
// 成分外の木を通る近道を使わず、両端から距離12の集合を分離する。
assert_eq!(
minimum_vertex_cut(
&input,
&vertices,
Coord::new(0, 0),
Coord::new(0, 2),
TERMINAL_RADIUS
),
[bridge]
);
}
#[test]
fn minimum_cut_may_remove_a_source_vertex() {
let mut input = Input {
d: Map2d::from_fn(|_| 5, Input::N),
};
let vertices: Vec<_> = (0..2 * TERMINAL_RADIUS + 4)
.map(|col| Coord::new(0, col))
.collect();
// 選択成分外にも木があるが、フローには含めない。
// sourceの右端を最安にし、その頂点自身を切れることも確認する。
let cheapest = Coord::new(0, TERMINAL_RADIUS);
input.d[cheapest] = 9;
assert_eq!(
minimum_vertex_cut(
&input,
&vertices,
vertices[0],
*vertices.last().unwrap(),
TERMINAL_RADIUS
),
[cheapest]
);
}
#[test]
fn removing_batches_restores_whole_cuts_and_preserves_best_snapshot() {
let mut state = State::new();
for row in 0..8 {
state.add((0..row + 1).map(|col| Coord::new(row, col)).collect());
}
let best = state.output();
let mut expected_rng = Random::new(42);
let expected_count =
REMOVE_LOWER + expected_rng.index(REMOVE_UPPER - REMOVE_LOWER + 1);
state.remove_random_batches(&mut Random::new(42));
assert_eq!(state.batches.len(), 8 - expected_count);
assert_eq!(best.cuts.len(), 36);
assert_eq!(state.count, state.cut.iter().filter(|&&cut| cut).count());
assert_eq!(state.count, state.output().cuts.len());
for row in 0..8 {
let retained = state.batches.iter().any(|batch| batch[0].row() == row);
for col in 0..row + 1 {
assert_eq!(state.cut[Coord::new(row, col)], retained);
}
}
}
#[test]
fn removing_more_batches_than_exist_restores_all_trees() {
let mut state = State::new();
state.add(vec![Coord::new(0, 0), Coord::new(0, 1)]);
state.remove_random_batches(&mut Random::new(42));
assert!(state.batches.is_empty());
assert_eq!(state.count, 0);
assert!(!state.cut.iter().any(|&value| value));
state.remove_random_batches(&mut Random::new(42));
assert_eq!(state.count, 0);
}
#[test]
fn minimum_cut_cost_matches_exhaustive_search() {
let mut input = Input {
d: Map2d::with_default(Input::N),
};
let vertices: Vec<_> = (0..12).map(|i| Coord::new(i / 4, i % 4)).collect();
for (i, &p) in vertices.iter().enumerate() {
input.d[p] = 5 + (i % 5) as u8;
}
let p1 = vertices[0];
let p2 = vertices[11];
let radius = 1;
let mut optimum = i64::MAX;
for mask in 0..(1 << vertices.len()) {
let cost: i64 = vertices
.iter()
.enumerate()
.filter(|(i, _)| mask & (1 << i) != 0)
.map(|(_, &p)| VERTEX_BASE_COST - i64::from(input.d[p]))
.sum();
if cost >= optimum {
continue;
}
let mut seen = vec![false; vertices.len()];
let mut stack = Vec::new();
for (i, &p) in vertices.iter().enumerate() {
if mask & (1 << i) == 0 && p.dist(&p1) <= radius {
seen[i] = true;
stack.push(i);
}
}
while let Some(i) = stack.pop() {
for (j, &p) in vertices.iter().enumerate() {
if mask & (1 << j) == 0 && !seen[j] && vertices[i].dist(&p) == 1 {
seen[j] = true;
stack.push(j);
}
}
}
if !vertices
.iter()
.enumerate()
.any(|(i, &p)| seen[i] && p.dist(&p2) <= radius)
{
optimum = cost;
}
}
let actual = minimum_vertex_cut(&input, &vertices, p1, p2, radius);
let cost: i64 = actual
.iter()
.map(|&p| VERTEX_BASE_COST - i64::from(input.d[p]))
.sum();
assert_eq!(cost, optimum);
}
}
#[cfg(test)]
fn components_reference(input: &Input, cut: &Map2d<bool>) -> Vec<Vec<Coord>> {
let mut seen = Map2d::with_default(Input::N);
let mut result = Vec::new();
for row in 0..Input::N {
for col in 0..Input::N {
let start = Coord::new(row, col);
if seen[start] || cut[start] || input.d[start] < D_THRESHOLD {
continue;
}
seen[start] = true;
let mut vertices = vec![start];
let mut head = 0;
while head < vertices.len() {
let p = vertices[head];
head += 1;
for diff in ADJACENTS {
let next = p + diff;
if next.in_map(Input::N)
&& !seen[next]
&& !cut[next]
&& input.d[next] >= D_THRESHOLD
{
seen[next] = true;
vertices.push(next);
}
}
}
result.push(vertices);
}
}
result
}
}
pub use vertex_cut::solve;
}
fn main() {
let timer = library::diagnostics::Timer::default();
let input = problem::Input::read();
let output = solver::solve(&input, &timer);
print!("{output}");
}