#1169446 (Rust) No.5007 Steiner Space Travel

• 問題: No.5007 Steiner Space Travel
• ユーザー: assy1028
• 提出日時: 2026-06-20 16:11:08
• 言語: Rust (1.94.0 + proconio + num + itertools)
• 結果: AC
• 実行時間: 703 ms / 1,000 ms
• コード長: 25,297 bytes
• コンパイル時間: 11,369 ms
• コンパイル使用メモリ: 210,204 KB
• 実行使用メモリ: 6,400 KB
• スコア: 8,962,864
• 最終ジャッジ日時: 2026-06-20 16:11:44
• 合計ジャッジ時間: 26,401 ms
• ジャッジサーバーID （参考情報）: judge1_1 / judge2_1

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

ソースコード

use std::io::{self, Read};
use std::time::{Duration, Instant};

const ALPHA: i64 = 5;
const COORD_MIN: i32 = 0;
const COORD_MAX: i32 = 1000;
const INF: i64 = 1_i64 << 60;

#[derive(Clone, Copy, Debug, PartialEq, Eq)]
struct Point {
    x: i32,
    y: i32,
}

#[derive(Clone)]
struct StationGraph {
    dist: Vec<Vec<i64>>,
    next: Vec<Vec<usize>>,
}

#[derive(Clone, Copy)]
struct PairPath {
    cost: i64,
    via_station: bool,
    s: usize,
    t: usize,
}

#[derive(Clone)]
struct Candidate {
    stations: Vec<Point>,
    order: Vec<usize>,
    cost: i64,
}

struct XorShift {
    state: u64,
}

impl XorShift {
    fn new(seed: u64) -> Self {
        Self { state: seed.max(1) }
    }

    fn next_u64(&mut self) -> u64 {
        let mut x = self.state;
        x ^= x << 7;
        x ^= x >> 9;
        self.state = x;
        x
    }

    fn gen_usize(&mut self, upper: usize) -> usize {
        if upper == 0 {
            0
        } else {
            (self.next_u64() as usize) % upper
        }
    }

    fn gen_f64(&mut self) -> f64 {
        const DEN: f64 = (1_u64 << 53) as f64;
        ((self.next_u64() >> 11) as f64) / DEN
    }

    fn gen_i32(&mut self, lo: i32, hi: i32) -> i32 {
        lo + self.gen_usize((hi - lo + 1) as usize) as i32
    }
}

fn dist2(a: Point, b: Point) -> i64 {
    let dx = (a.x - b.x) as i64;
    let dy = (a.y - b.y) as i64;
    dx * dx + dy * dy
}

fn clamp_coord(v: f64) -> i32 {
    (v.round() as i32).clamp(COORD_MIN, COORD_MAX)
}

fn mean_point(points: &[Point], ids: &[usize]) -> Point {
    if ids.is_empty() {
        return Point { x: 500, y: 500 };
    }
    let mut sx = 0.0;
    let mut sy = 0.0;
    for &id in ids {
        sx += points[id].x as f64;
        sy += points[id].y as f64;
    }
    let len = ids.len() as f64;
    Point {
        x: clamp_coord(sx / len),
        y: clamp_coord(sy / len),
    }
}

fn build_station_graph(stations: &[Point]) -> StationGraph {
    let m = stations.len();
    let mut dist = vec![vec![0; m]; m];
    let mut next = vec![vec![0; m]; m];
    for i in 0..m {
        for j in 0..m {
            dist[i][j] = dist2(stations[i], stations[j]);
            next[i][j] = j;
        }
    }
    for k in 0..m {
        for i in 0..m {
            for j in 0..m {
                let nd = dist[i][k] + dist[k][j];
                if nd < dist[i][j] {
                    dist[i][j] = nd;
                    next[i][j] = next[i][k];
                }
            }
        }
    }
    StationGraph { dist, next }
}

fn planet_pair_path(
    points: &[Point],
    stations: &[Point],
    sg: &StationGraph,
    a: usize,
    b: usize,
) -> PairPath {
    let mut best = PairPath {
        cost: ALPHA * ALPHA * dist2(points[a], points[b]),
        via_station: false,
        s: 0,
        t: 0,
    };
    for s in 0..stations.len() {
        let from_a = ALPHA * dist2(points[a], stations[s]);
        for t in 0..stations.len() {
            let cost = from_a + sg.dist[s][t] + ALPHA * dist2(stations[t], points[b]);
            if cost < best.cost {
                best = PairPath {
                    cost,
                    via_station: true,
                    s,
                    t,
                };
            }
        }
    }
    best
}

fn build_planet_dist(points: &[Point], stations: &[Point]) -> Vec<Vec<i64>> {
    let n = points.len();
    let sg = build_station_graph(stations);
    let mut d = vec![vec![0; n]; n];
    for i in 0..n {
        for j in (i + 1)..n {
            let c = planet_pair_path(points, stations, &sg, i, j).cost;
            d[i][j] = c;
            d[j][i] = c;
        }
    }
    d
}

fn route_cost(order: &[usize], dist: &[Vec<i64>]) -> i64 {
    let n = order.len();
    if n == 0 {
        return 0;
    }
    let mut cost = 0;
    for i in 0..n {
        cost += dist[order[i]][order[(i + 1) % n]];
    }
    cost
}

fn nearest_neighbor(dist: &[Vec<i64>], first: Option<usize>) -> Vec<usize> {
    let n = dist.len();
    if n == 0 {
        return Vec::new();
    }
    let mut used = vec![false; n];
    let mut order = vec![0];
    used[0] = true;
    if let Some(f) = first {
        if f != 0 && f < n {
            order.push(f);
            used[f] = true;
        }
    }
    while order.len() < n {
        let cur = *order.last().unwrap();
        let mut best = None;
        for v in 0..n {
            if !used[v] && best.map_or(true, |b| dist[cur][v] < dist[cur][b]) {
                best = Some(v);
            }
        }
        let v = best.unwrap();
        used[v] = true;
        order.push(v);
    }
    order
}

fn farthest_insertion(dist: &[Vec<i64>]) -> Vec<usize> {
    let n = dist.len();
    if n <= 1 {
        return (0..n).collect();
    }
    let mut used = vec![false; n];
    let mut route = vec![0];
    used[0] = true;
    while route.len() < n {
        let mut pick = 0;
        let mut pick_score = -1_i64;
        for v in 0..n {
            if used[v] {
                continue;
            }
            let nearest = route.iter().map(|&r| dist[v][r]).min().unwrap();
            if nearest > pick_score {
                pick_score = nearest;
                pick = v;
            }
        }
        insert_at_best_position(&mut route, pick, dist);
        used[pick] = true;
    }
    route
}

fn cheapest_insertion(dist: &[Vec<i64>]) -> Vec<usize> {
    let n = dist.len();
    if n <= 1 {
        return (0..n).collect();
    }
    let mut used = vec![false; n];
    let mut route = vec![0];
    used[0] = true;
    while route.len() < n {
        let mut best_v = 0;
        let mut best_pos = 0;
        let mut best_delta = INF;
        for v in 0..n {
            if used[v] {
                continue;
            }
            let len = route.len();
            for pos in 0..len {
                let a = route[pos];
                let b = route[(pos + 1) % len];
                let delta = dist[a][v] + dist[v][b] - dist[a][b];
                if delta < best_delta {
                    best_delta = delta;
                    best_v = v;
                    best_pos = pos + 1;
                }
            }
        }
        route.insert(best_pos, best_v);
        used[best_v] = true;
    }
    route
}

fn insert_at_best_position(route: &mut Vec<usize>, v: usize, dist: &[Vec<i64>]) {
    let len = route.len();
    let mut best_pos = 1;
    let mut best_delta = INF;
    for pos in 0..len {
        let a = route[pos];
        let b = route[(pos + 1) % len];
        let delta = dist[a][v] + dist[v][b] - dist[a][b];
        if delta < best_delta {
            best_delta = delta;
            best_pos = pos + 1;
        }
    }
    route.insert(best_pos, v);
}

fn two_opt(order: &mut [usize], dist: &[Vec<i64>]) -> bool {
    let n = order.len();
    if n <= 3 {
        return false;
    }
    let mut any = false;
    for _ in 0..25 {
        let mut improved = false;
        for i in 0..(n - 1) {
            let a = order[i];
            let b = order[(i + 1) % n];
            for j in (i + 2)..n {
                if i == 0 && j == n - 1 {
                    continue;
                }
                let c = order[j];
                let d = order[(j + 1) % n];
                let delta = dist[a][c] + dist[b][d] - dist[a][b] - dist[c][d];
                if delta < 0 {
                    order[(i + 1)..=j].reverse();
                    improved = true;
                    any = true;
                }
            }
        }
        if !improved {
            break;
        }
    }
    any
}

fn relocate(order: &mut Vec<usize>, dist: &[Vec<i64>]) -> bool {
    let n = order.len();
    if n <= 3 {
        return false;
    }
    for _ in 0..20 {
        let mut improved = false;
        'outer: for i in 1..n {
            let prev = order[i - 1];
            let x = order[i];
            let next = order[(i + 1) % n];
            let remove_delta = dist[prev][next] - dist[prev][x] - dist[x][next];
            for j in 0..n {
                if j == i || (j + 1) % n == i {
                    continue;
                }
                let a = order[j];
                let b = order[(j + 1) % n];
                let delta = remove_delta + dist[a][x] + dist[x][b] - dist[a][b];
                if delta < 0 {
                    let moved = order.remove(i);
                    let insert_pos = if j < i { j + 1 } else { j };
                    order.insert(insert_pos, moved);
                    improved = true;
                    break 'outer;
                }
            }
        }
        if !improved {
            return false;
        }
    }
    true
}

fn improve_order(order: &mut Vec<usize>, dist: &[Vec<i64>]) {
    for _ in 0..4 {
        let a = two_opt(order, dist);
        let b = relocate(order, dist);
        if !a && !b {
            break;
        }
    }
}

fn solve_tsp(dist: &[Vec<i64>], rng: &mut XorShift, extra_random: usize) -> Vec<usize> {
    let n = dist.len();
    if n <= 1 {
        return (0..n).collect();
    }
    let mut best = farthest_insertion(dist);
    improve_order(&mut best, dist);
    let mut best_cost = route_cost(&best, dist);

    let mut candidates = Vec::new();
    let mut by_start: Vec<usize> = (1..n).collect();
    by_start.sort_by_key(|&v| dist[0][v]);
    candidates.extend(by_start.iter().take(2).copied());
    candidates.extend(by_start.iter().rev().take(2).copied());
    candidates.push(usize::MAX);
    for _ in 0..extra_random {
        candidates.push(1 + rng.gen_usize(n - 1));
    }

    for first in candidates {
        let mut order = if first == usize::MAX {
            cheapest_insertion(dist)
        } else {
            nearest_neighbor(dist, Some(first))
        };
        improve_order(&mut order, dist);
        let cost = route_cost(&order, dist);
        if cost < best_cost {
            best_cost = cost;
            best = order;
        }
    }
    best
}

fn evaluate_candidate(points: &[Point], stations: Vec<Point>, rng: &mut XorShift) -> Candidate {
    let dist = build_planet_dist(points, &stations);
    let order = solve_tsp(&dist, rng, 3);
    let cost = route_cost(&order, &dist);
    Candidate {
        stations,
        order,
        cost,
    }
}

fn fixed_order_cost(points: &[Point], stations: &[Point], order: &[usize]) -> i64 {
    if order.is_empty() {
        return 0;
    }
    let sg = build_station_graph(stations);
    let mut cost = 0;
    for i in 0..order.len() {
        cost += planet_pair_path(
            points,
            stations,
            &sg,
            order[i],
            order[(i + 1) % order.len()],
        )
        .cost;
    }
    cost
}

fn coordinate_descent(points: &[Point], cand: &mut Candidate, deadline: Instant) {
    let dirs = [
        (1, 0),
        (-1, 0),
        (0, 1),
        (0, -1),
        (1, 1),
        (1, -1),
        (-1, 1),
        (-1, -1),
    ];
    let steps = [96, 64, 40, 24, 16, 10, 6, 3, 1];
    let mut cur_cost = fixed_order_cost(points, &cand.stations, &cand.order);
    for &step in &steps {
        loop {
            if Instant::now() >= deadline {
                cand.cost = cur_cost;
                return;
            }
            let mut improved = false;
            for s in 0..cand.stations.len() {
                let orig = cand.stations[s];
                let mut best_point = orig;
                let mut best_cost = cur_cost;
                for &(dx, dy) in &dirs {
                    let p = Point {
                        x: (orig.x + dx * step).clamp(COORD_MIN, COORD_MAX),
                        y: (orig.y + dy * step).clamp(COORD_MIN, COORD_MAX),
                    };
                    if p == orig {
                        continue;
                    }
                    cand.stations[s] = p;
                    let cost = fixed_order_cost(points, &cand.stations, &cand.order);
                    if cost < best_cost {
                        best_cost = cost;
                        best_point = p;
                    }
                }
                cand.stations[s] = best_point;
                if best_cost < cur_cost {
                    cur_cost = best_cost;
                    improved = true;
                }
                if Instant::now() >= deadline {
                    cand.cost = cur_cost;
                    return;
                }
            }
            if !improved {
                break;
            }
        }
    }
    cand.cost = cur_cost;
}

fn refine_candidate(
    points: &[Point],
    mut cand: Candidate,
    rng: &mut XorShift,
    deadline: Instant,
) -> Candidate {
    coordinate_descent(points, &mut cand, deadline);
    if Instant::now() < deadline {
        let dist = build_planet_dist(points, &cand.stations);
        improve_order(&mut cand.order, &dist);
        cand.cost = route_cost(&cand.order, &dist);
    }
    if Instant::now() < deadline {
        coordinate_descent(points, &mut cand, deadline);
    }
    if Instant::now() < deadline {
        let dist = build_planet_dist(points, &cand.stations);
        let order = solve_tsp(&dist, rng, 6);
        let cost = route_cost(&order, &dist);
        if cost < cand.cost {
            cand.order = order;
            cand.cost = cost;
        }
    }
    cand
}

fn lloyd(points: &[Point], mut centers: Vec<(f64, f64)>, iters: usize) -> Vec<Point> {
    let k = centers.len();
    if k == 0 {
        return Vec::new();
    }
    let n = points.len();
    let mut assign = vec![0; n];
    for _ in 0..iters {
        for (i, &p) in points.iter().enumerate() {
            let mut best = 0;
            let mut best_d = f64::INFINITY;
            for (c, &(x, y)) in centers.iter().enumerate() {
                let dx = p.x as f64 - x;
                let dy = p.y as f64 - y;
                let d = dx * dx + dy * dy;
                if d < best_d {
                    best_d = d;
                    best = c;
                }
            }
            assign[i] = best;
        }

        let mut sx = vec![0.0; k];
        let mut sy = vec![0.0; k];
        let mut cnt = vec![0_usize; k];
        for (i, &a) in assign.iter().enumerate() {
            sx[a] += points[i].x as f64;
            sy[a] += points[i].y as f64;
            cnt[a] += 1;
        }

        for c in 0..k {
            if cnt[c] > 0 {
                centers[c] = (sx[c] / cnt[c] as f64, sy[c] / cnt[c] as f64);
            } else {
                let mut farthest = 0;
                let mut best_d = -1.0;
                for (i, &p) in points.iter().enumerate() {
                    let d = centers
                        .iter()
                        .map(|&(x, y)| {
                            let dx = p.x as f64 - x;
                            let dy = p.y as f64 - y;
                            dx * dx + dy * dy
                        })
                        .fold(f64::INFINITY, f64::min);
                    if d > best_d {
                        best_d = d;
                        farthest = i;
                    }
                }
                centers[c] = (points[farthest].x as f64, points[farthest].y as f64);
            }
        }
    }
    centers
        .into_iter()
        .map(|(x, y)| Point {
            x: clamp_coord(x),
            y: clamp_coord(y),
        })
        .collect()
}

fn kmeans_plus_plus(points: &[Point], k: usize, rng: &mut XorShift) -> Vec<Point> {
    if k == 0 {
        return Vec::new();
    }
    let n = points.len();
    let mut centers = Vec::with_capacity(k);
    let first = rng.gen_usize(n);
    centers.push((points[first].x as f64, points[first].y as f64));
    while centers.len() < k {
        let mut weights = vec![0.0; n];
        let mut sum = 0.0;
        for (i, &p) in points.iter().enumerate() {
            let d = centers
                .iter()
                .map(|&(x, y)| {
                    let dx = p.x as f64 - x;
                    let dy = p.y as f64 - y;
                    dx * dx + dy * dy
                })
                .fold(f64::INFINITY, f64::min);
            weights[i] = d;
            sum += d;
        }
        let pick = if sum <= 0.0 {
            rng.gen_usize(n)
        } else {
            let mut r = rng.gen_f64() * sum;
            let mut id = n - 1;
            for (i, &w) in weights.iter().enumerate() {
                r -= w;
                if r <= 0.0 {
                    id = i;
                    break;
                }
            }
            id
        };
        centers.push((points[pick].x as f64, points[pick].y as f64));
    }
    lloyd(points, centers, 20)
}

fn farthest_kmeans(points: &[Point], k: usize) -> Vec<Point> {
    if k == 0 {
        return Vec::new();
    }
    let n = points.len();
    let mut sx = 0.0;
    let mut sy = 0.0;
    for &p in points {
        sx += p.x as f64;
        sy += p.y as f64;
    }
    let centroid = Point {
        x: clamp_coord(sx / n as f64),
        y: clamp_coord(sy / n as f64),
    };
    let mut first = 0;
    let mut best = INF;
    for i in 0..n {
        let d = dist2(points[i], centroid);
        if d < best {
            best = d;
            first = i;
        }
    }
    let mut ids = vec![first];
    while ids.len() < k {
        let mut pick = 0;
        let mut best_d = -1_i64;
        for i in 0..n {
            let d = ids
                .iter()
                .map(|&c| dist2(points[i], points[c]))
                .min()
                .unwrap();
            if d > best_d {
                best_d = d;
                pick = i;
            }
        }
        ids.push(pick);
    }
    let centers = ids
        .iter()
        .map(|&i| (points[i].x as f64, points[i].y as f64))
        .collect();
    lloyd(points, centers, 20)
}

fn chunk_stations(points: &[Point], k: usize, key: usize) -> Vec<Point> {
    let mut ids: Vec<usize> = (0..points.len()).collect();
    ids.sort_by_key(|&i| match key {
        0 => points[i].x,
        1 => points[i].y,
        2 => points[i].x + points[i].y,
        _ => points[i].x - points[i].y,
    });
    let mut stations = Vec::with_capacity(k);
    for c in 0..k {
        let l = c * ids.len() / k;
        let r = ((c + 1) * ids.len() / k).max(l + 1).min(ids.len());
        stations.push(mean_point(points, &ids[l..r]));
    }
    let centers = stations.iter().map(|p| (p.x as f64, p.y as f64)).collect();
    lloyd(points, centers, 12)
}

fn grid_stations(points: &[Point], k: usize, cols: usize) -> Vec<Point> {
    let rows = (k + cols - 1) / cols;
    let min_x = points.iter().map(|p| p.x).min().unwrap_or(0) as f64;
    let max_x = points.iter().map(|p| p.x).max().unwrap_or(1000) as f64;
    let min_y = points.iter().map(|p| p.y).min().unwrap_or(0) as f64;
    let max_y = points.iter().map(|p| p.y).max().unwrap_or(1000) as f64;
    let mut stations = Vec::with_capacity(k);
    for r in 0..rows {
        for c in 0..cols {
            if stations.len() == k {
                break;
            }
            let x = min_x + (max_x - min_x) * (c as f64 + 0.5) / cols as f64;
            let y = min_y + (max_y - min_y) * (r as f64 + 0.5) / rows as f64;
            stations.push(Point {
                x: clamp_coord(x),
                y: clamp_coord(y),
            });
        }
    }
    stations
}

fn random_perturb(stations: &[Point], rng: &mut XorShift, radius: i32) -> Vec<Point> {
    let mut res = stations.to_vec();
    for p in &mut res {
        p.x = (p.x + rng.gen_i32(-radius, radius)).clamp(COORD_MIN, COORD_MAX);
        p.y = (p.y + rng.gen_i32(-radius, radius)).clamp(COORD_MIN, COORD_MAX);
    }
    res
}

fn input_seed(points: &[Point], m: usize) -> u64 {
    let mut h = 0x9e37_79b9_7f4a_7c15_u64 ^ m as u64;
    for &p in points {
        h ^= (p.x as u64).wrapping_mul(0xbf58_476d_1ce4_e5b9);
        h = h.rotate_left(13);
        h ^= (p.y as u64).wrapping_mul(0x94d0_49bb_1331_11eb);
        h = h.rotate_left(17);
    }
    h
}

fn append_segment(
    points: &[Point],
    stations: &[Point],
    sg: &StationGraph,
    a: usize,
    b: usize,
    out: &mut Vec<(i32, usize)>,
) {
    let path = planet_pair_path(points, stations, sg, a, b);
    if path.via_station {
        out.push((2, path.s + 1));
        let mut cur = path.s;
        while cur != path.t {
            cur = sg.next[cur][path.t];
            out.push((2, cur + 1));
        }
    }
    out.push((1, b + 1));
}

fn build_visit_sequence(
    points: &[Point],
    stations: &[Point],
    order: &[usize],
) -> Vec<(i32, usize)> {
    let mut out = Vec::new();
    if order.is_empty() {
        return out;
    }
    let sg = build_station_graph(stations);
    out.push((1, 1));
    for i in 0..order.len() {
        append_segment(
            points,
            stations,
            &sg,
            order[i],
            order[(i + 1) % order.len()],
            &mut out,
        );
    }
    out
}

fn solve(points: &[Point], m: usize) -> Candidate {
    let start = Instant::now();
    let generation_deadline = start + Duration::from_millis(430);
    let refine_deadline = start + Duration::from_millis(700);
    let mut rng = XorShift::new(input_seed(points, m));

    let mut best: Option<Candidate> = None;
    let try_candidate = |stations: Vec<Point>, rng: &mut XorShift, best: &mut Option<Candidate>| {
        let cand = evaluate_candidate(points, stations, rng);
        if best.as_ref().map_or(true, |b| cand.cost < b.cost) {
            *best = Some(cand);
        }
    };

    if m == 0 {
        let stations = Vec::new();
        let dist = build_planet_dist(points, &stations);
        let order = solve_tsp(&dist, &mut rng, 3);
        let cost = route_cost(&order, &dist);
        return Candidate {
            stations,
            order,
            cost,
        };
    }

    try_candidate(farthest_kmeans(points, m), &mut rng, &mut best);
    for key in 0..4 {
        try_candidate(chunk_stations(points, m, key), &mut rng, &mut best);
    }
    for cols in 1..=m {
        try_candidate(grid_stations(points, m, cols), &mut rng, &mut best);
    }

    let mut iter = 0;
    while Instant::now() < generation_deadline && iter < 120 {
        let stations = if iter % 5 == 0 {
            let base = best
                .as_ref()
                .map(|b| b.stations.clone())
                .unwrap_or_else(|| kmeans_plus_plus(points, m, &mut rng));
            random_perturb(&base, &mut rng, 70)
        } else {
            kmeans_plus_plus(points, m, &mut rng)
        };
        try_candidate(stations, &mut rng, &mut best);
        iter += 1;
    }

    let mut cand = best.unwrap();
    cand = refine_candidate(points, cand, &mut rng, refine_deadline);

    while Instant::now() < refine_deadline {
        let radius = 30.max(120 - (start.elapsed().as_millis() as i32 / 6));
        let stations = random_perturb(&cand.stations, &mut rng, radius);
        let mut next = evaluate_candidate(points, stations, &mut rng);
        if Instant::now() < refine_deadline {
            coordinate_descent(points, &mut next, refine_deadline);
            let dist = build_planet_dist(points, &next.stations);
            improve_order(&mut next.order, &dist);
            next.cost = route_cost(&next.order, &dist);
        }
        if next.cost < cand.cost {
            cand = next;
        }
    }

    cand
}

fn parse_input() -> (usize, usize, Vec<Point>) {
    let mut s = String::new();
    io::stdin().read_to_string(&mut s).unwrap();
    let mut it = s.split_whitespace();
    let n: usize = it.next().unwrap().parse().unwrap();
    let m: usize = it.next().unwrap().parse().unwrap();
    let mut points = Vec::with_capacity(n);
    for _ in 0..n {
        let x: i32 = it.next().unwrap().parse().unwrap();
        let y: i32 = it.next().unwrap().parse().unwrap();
        points.push(Point { x, y });
    }
    (n, m, points)
}

fn main() {
    let (_n, m, points) = parse_input();
    let cand = solve(&points, m);
    let route = build_visit_sequence(&points, &cand.stations, &cand.order);

    let mut out = String::new();
    for p in &cand.stations {
        out.push_str(&format!("{} {}\n", p.x, p.y));
    }
    out.push_str(&format!("{}\n", route.len()));
    for (t, r) in route {
        out.push_str(&format!("{} {}\n", t, r));
    }
    print!("{out}");
}

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

    #[test]
    fn one_station_can_reduce_two_planet_edge() {
        let points = vec![Point { x: 0, y: 0 }, Point { x: 200, y: 200 }];
        let stations = vec![Point { x: 100, y: 100 }];
        let sg = build_station_graph(&stations);
        let path = planet_pair_path(&points, &stations, &sg, 0, 1);
        assert!(path.via_station);
        assert_eq!(path.cost, 200_000);
    }

    #[test]
    fn output_sequence_starts_and_ends_at_planet_one() {
        let points = vec![
            Point { x: 0, y: 0 },
            Point { x: 200, y: 200 },
            Point { x: 0, y: 200 },
        ];
        let stations = vec![Point { x: 100, y: 100 }];
        let order = vec![0, 1, 2];
        let seq = build_visit_sequence(&points, &stations, &order);
        assert_eq!(seq.first(), Some(&(1, 1)));
        assert_eq!(seq.last(), Some(&(1, 1)));
        for i in 1..=3 {
            assert!(seq.contains(&(1, i)));
        }
    }
}