#![allow(non_snake_case)] use std::collections::BinaryHeap; fn main() { get_time(); let input = read_input(); let out = solve(&input); eprintln!("Time = {:.3}", get_time()); eprintln!("Score = {}", compute_score(&input, &out)); write_output(&input, &out); // solution::write_profile(); } const INF: i64 = 1000000000000000; fn solve(input: &Input) -> Output { let mut rng = XorShift::new(4823902); let mut ps = input.ps.clone(); for _ in 0..input.M { ps.push((rng.next(801) as i32 + 100, rng.next(801) as i32 + 100)); } let mut out = Output { ps, route: vec![] }; let mut best = 0; let mut rep = 0; while get_time() < 0.9 { let last = get_time() > 0.85; rep += 1; let mut ps = out.ps.clone(); if !last { let p = rng.next(input.M); ps[input.N + p] = (rng.next(801) as i32 + 100, rng.next(801) as i32 + 100); } let mut g = mat![0; input.N; input.N]; let mut inter = mat![!0; input.N; input.N]; for i in 0..input.N { for j in 0..i { g[i][j] = dist2(ps[i], ps[j]) * 5 * 5; g[j][i] = g[i][j]; } } let mut prevs = vec![vec![]; input.N]; for i in 0..input.N { let mut dist = vec![INF; input.M]; let mut prev = vec![!0; input.M]; let mut que = BinaryHeap::new(); for j in 0..input.M { dist[j] = dist2(ps[i], ps[input.N + j]) * 5; que.push((-dist[j], j)); } while let Some((d, u)) = que.pop() { let d = -d; if d != dist[u] { continue; } for v in 0..input.M { if dist[v].setmin(d + dist2(ps[input.N + u], ps[input.N + v])) { prev[v] = u; que.push((-dist[v], v)); } } } for j in 0..input.N { for t in 0..input.M { let d = dist[t] + dist2(ps[input.N + t], ps[j]) * 5; if g[i][j].setmin(d) { inter[i][j] = t; } } } prevs[i] = prev; } let mut route = if out.route.is_empty() { tsp::greedy(&g) } else { out.route.iter().filter(|&&v| v < input.N).cloned().collect() }; route = tsp::solve(&g, &route, if last { 0.95 } else { get_time() }, &mut rng); let mut route2 = vec![route[0]]; for i in 1..route.len() { if inter[route[i - 1]][route[i]] != !0 { let mut v = inter[route[i - 1]][route[i]]; let mut tmp = vec![]; while v != !0 { tmp.push(input.N + v); v = prevs[route[i - 1]][v]; } tmp.reverse(); route2.extend(tmp); } route2.push(route[i]); } let crt = Output { ps: ps.clone(), route: route2.clone() }; if best.setmax(compute_score(input, &crt)) { out = crt; eprintln!("{:.3}: {}", get_time(), best); } for xy in 0..2 { let mut c1 = vec![0; input.M]; let mut c2 = mat![0; input.M; input.M]; for i in 1..route2.len() { if route2[i - 1] < input.N && route2[i] >= input.N { let a = route2[i] - input.N; let c = if xy == 0 { ps[route2[i - 1]].0 } else { ps[route2[i - 1]].1 }; c2[a][a] += 5; c1[a] -= 10 * c; } else if route2[i - 1] >= input.N && route2[i] < input.N { let a = route2[i - 1] - input.N; let c = if xy == 0 { ps[route2[i]].0 } else { ps[route2[i]].1 }; c2[a][a] += 5; c1[a] -= 10 * c; } else if route2[i - 1] >= input.N && route2[i] >= input.N { let a = route2[i - 1] - input.N; let b = route2[i] - input.N; if a != b { c1[a] += 1; c1[b] += 1; c2[a][b] -= 2; c2[b][a] -= 2; } } } let mut A = mat![0.0; input.M; input.M]; let mut b = vec![0.0; input.M]; for i in 0..input.M { for j in 0..input.M { A[i][j] = c2[i][j] as f64; } A[i][i] *= 2.0; b[i] = -c1[i] as f64; } let A = Mat(A); let decomp = A.gauss(); let x = decomp.solve(b); if x.len() > 0 { let x = x[0].clone(); for i in 0..input.M { if xy == 0 { ps[input.N + i].0 = x[i].round() as i32; } else { ps[input.N + i].1 = x[i].round() as i32; } } } } let crt = Output { ps: ps.clone(), route: route2.clone() }; if best.setmax(compute_score(input, &crt)) { out = crt; eprintln!(" {:.3}: {}", get_time(), best); } } eprintln!("rep = {}", rep); out } // 入出力と得点計算 struct Output { ps: Vec<(i32, i32)>, route: Vec, } struct Input { N: usize, M: usize, ps: Vec<(i32, i32)>, } fn read_input() -> Input { let (N, M) = read!(usize, usize); let ps = read!(i32, i32; N); Input { N, M, ps } } fn write_output(input: &Input, out: &Output) { for &(x, y) in &out.ps[input.N..] { println!("{} {}", x, y); } println!("{}", out.route.len()); for &p in &out.route { if p < input.N { println!("1 {}", p + 1); } else { println!("2 {}", p - input.N + 1); } } } fn dist2((x1, y1): (i32, i32), (x2, y2): (i32, i32)) -> i64 { ((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)) as i64 } fn compute_score(input: &Input, out: &Output) -> i64 { let mut cost = 0; for i in 1..out.route.len() { let mut tmp = dist2(out.ps[out.route[i - 1]], out.ps[out.route[i]]); if out.route[i - 1] < input.N { tmp *= 5; } if out.route[i] < input.N { tmp *= 5; } cost += tmp; } (1e9 / (1e3 + (cost as f64).sqrt())).round() as i64 } // ここからライブラリ pub trait SetMinMax { fn setmin(&mut self, v: Self) -> bool; fn setmax(&mut self, v: Self) -> bool; } impl SetMinMax for T where T: PartialOrd { fn setmin(&mut self, v: T) -> bool { *self > v && { *self = v; true } } fn setmax(&mut self, v: T) -> bool { *self < v && { *self = v; true } } } #[macro_export] macro_rules! mat { ($($e:expr),*) => { vec![$($e),*] }; ($($e:expr,)*) => { vec![$($e),*] }; ($e:expr; $d:expr) => { vec![$e; $d] }; ($e:expr; $d:expr $(; $ds:expr)+) => { vec![mat![$e $(; $ds)*]; $d] }; } pub fn get_time() -> f64 { static mut STIME: f64 = -1.0; let t = std::time::SystemTime::now().duration_since(std::time::UNIX_EPOCH).unwrap(); let ms = t.as_secs() as f64 + t.subsec_nanos() as f64 * 1e-9; unsafe { if STIME < 0.0 { STIME = ms; } // ローカル環境とジャッジ環境の実行速度差はget_timeで吸収しておくと便利 #[cfg(feature="local")] { (ms - STIME) * 2.0 } #[cfg(not(feature="local"))] { (ms - STIME) } } } mod tsp { use super::*; pub fn compute_cost(g: &Vec>, ps: &Vec) -> i64 { let mut tmp = 0; for i in 0..ps.len() - 1 { tmp += g[ps[i]][ps[i + 1]]; } tmp } pub fn greedy(g: &Vec>) -> Vec { let mut ps = vec![0]; let n = g.len(); let mut used = vec![false; n]; used[0] = true; for i in 0..n-1 { let mut to = !0; let mut cost = i64::max_value(); for j in 0..n { if !used[j] && cost.setmin(g[i][j]) { to = j; } } used[to] = true; ps.push(to); } ps.push(0); ps } // mv: (i, dir) pub fn apply_move(tour: &mut Vec, idx: &mut Vec, mv: &[(usize, usize)]) { let k = mv.len(); let mut ids: Vec<_> = (0..k).collect(); ids.sort_by_key(|&i| mv[i].0); let mut order = vec![0; k]; for i in 0..k { order[ids[i]] = i; } let mut tour2 = Vec::with_capacity(mv[ids[k - 1]].0 - mv[ids[0]].0); let mut i = ids[0]; let mut dir = 0; loop { let (j, rev) = if dir == mv[i].1 { ((i + 1) % k, 0) } else { ((i + k - 1) % k, 1) }; if mv[j].1 == rev { if order[j] == k - 1 { break; } else { i = ids[order[j] + 1]; dir = 0; tour2.extend_from_slice(&tour[mv[j].0 + 1..mv[i].0 + 1]); } } else { i = ids[order[j] - 1]; dir = 1; tour2.extend(tour[mv[i].0 + 1..mv[j].0 + 1].iter().rev().cloned()); } } assert_eq!(tour2.len(), mv[ids[k - 1]].0 - mv[ids[0]].0); tour[mv[ids[0]].0 + 1..mv[ids[k - 1]].0 + 1].copy_from_slice(&tour2); for i in mv[ids[0]].0 + 1..mv[ids[k - 1]].0 + 1 { idx[tour[i]] = i; } } pub const FEASIBLE3: [bool; 64] = [false, false, false, true, false, true, true, true, true, true, true, false, true, false, false, false, false, false, false, false, false, false, false, false, false, false, false, true, false, true, true, true, true, true, true, false, true, false, false, false, false, false, false, false, false, false, false, false, false, false, false, true, false, true, true, true, true, true, true, false, true, false, false, false]; pub fn solve(g: &Vec>, qs: &Vec, until: f64, rng: &mut XorShift) -> Vec { // solution::profile!(tsp); let n = g.len(); let mut f = vec![vec![]; n]; for i in 0..n { for j in 0..n { if i != j { f[i].push((g[i][j], j)); } } f[i].sort_by(|&(a, _), &(b, _)| a.partial_cmp(&b).unwrap()); } let mut ps = qs.clone(); let mut idx = vec![!0; n]; let (mut min, mut min_ps) = (compute_cost(&g, &qs), ps.clone()); let mut first = true; while first || get_time() < until { first = false; let mut cost = compute_cost(&g, &ps); for p in 0..n { idx[ps[p]] = p; } loop { let mut ok = false; for i in 0..n { for di in 0..2 { 'loop_ij: for &(ij, vj) in &f[ps[i + di]] { if g[ps[i]][ps[i + 1]] - ij <= 0 { break; } for dj in 0..2 { let j = if idx[vj] == 0 && dj == 0 { n - 1 } else { idx[vj] - 1 + dj }; let gain = g[ps[i]][ps[i + 1]] - ij + g[ps[j]][ps[j + 1]]; // 2-opt if di != dj && gain - g[ps[j + dj]][ps[i + 1 - di]] > 0 { cost -= gain - g[ps[j + dj]][ps[i + 1 - di]]; apply_move(&mut ps, &mut idx, &[(i, di), (j, dj)]); ok = true; break 'loop_ij; } // 3-opt for &(jk, vk) in &f[ps[j + dj]] { if gain - jk <= 0 { break; } for dk in 0..2 { let k = if idx[vk] == 0 && dk == 0 { n - 1 } else { idx[vk] - 1 + dk }; if i == k || j == k { continue; } let gain = gain - jk + g[ps[k]][ps[k + 1]]; if gain - g[ps[k + dk]][ps[i + 1 - di]] > 0 { let mask = if i < j { 1 << 5 } else { 0 } | if i < k { 1 << 4 } else { 0 } | if j < k { 1 << 3 } else { 0 } | di << 2 | dj << 1 | dk; if FEASIBLE3[mask] { cost -= gain - g[ps[k + dk]][ps[i + 1 - di]]; apply_move(&mut ps, &mut idx, &[(i, di), (j, dj), (k, dk)]); ok = true; break 'loop_ij; } } } } } } } } if !ok { break; } } if min.setmin(cost) { min_ps = ps; } ps = min_ps.clone(); if n <= 4 { break; } loop { if rng.next(2) == 0 { // double bridge let mut is: Vec<_> = (0..4).map(|_| rng.next(n)).collect(); is.sort(); if is[0] == is[1] || is[1] == is[2] || is[2] == is[3] { continue; } ps = ps[0..is[0]+1].iter() .chain(ps[is[2]+1..is[3]+1].iter()) .chain(ps[is[1]+1..is[2]+1].iter()) .chain(ps[is[0]+1..is[1]+1].iter()) .chain(ps[is[3]+1..].iter()) .cloned().collect(); } else { for _ in 0..6 { loop { let i = 1 + rng.next(n - 1) as usize; let j = 1 + rng.next(n - 1) as usize; if i < j && j - i < n - 2 { ps = ps[0..i].iter().chain(ps[i..j+1].iter().rev()).chain(ps[j+1..].iter()).cloned().collect(); break; } } } } break; } } min_ps } } use std::ops::*; #[derive(Clone, Debug, Eq, Hash, Ord, PartialEq, PartialOrd)] pub struct Mat(pub Vec>); impl Mat where T: Clone + Default { pub fn row(a: Vec) -> Mat { Mat(vec![a]) } pub fn col(a: Vec) -> Mat { Mat(a.into_iter().map(|v| vec![v]).collect()) } pub fn to_vec(&self) -> Vec { if self.len() == 1 { self[0].clone() } else if self[0].len() == 1 { self.iter().map(|r| r[0].clone()).collect() } else { panic!("Cannot transform Mat into Vec. row: {}, col: {}", self.len(), self[0].len()); } } pub fn transpose(&self) -> Mat { let mut v = vec![Vec::with_capacity(self.len()); self[0].len()]; for r in self.iter() { for (i, x) in r.iter().enumerate() { v[i].push(x.clone()); } } Mat(v) } pub fn pow(&self, mut b: u64) -> Mat where T: AddAssign + PartialEq + From, for<'b> &'b T: Mul { let n = self.len(); assert_eq!(n, self[0].len()); let mut res: Mat = Mat(mat![0.into(); n; n]); for i in 0..n { res[i][i] = 1.into() } let mut a = self.clone(); while b > 0 { if (b & 1) > 0 { res = &res * &a; } a = &a * &a; b >>= 1; } res } /// O(nmr), r: rank. pub fn gauss(self) -> Decomposed where T: Clone + Default + PartialOrd + From + SubAssign, for<'b> &'b T: Neg + Mul + Div { let mut a = self; let (n, m) = (a.len(), a[0].len()); let (zero, one) = (T::default(), T::from(1)); let abs = |a: &T| { if a < &zero { a.neg() } else { a.clone() }}; let mut pivots = vec![n; m]; let mut r = 0; for c in 0..m { let mut t = abs(&a[r][c]); let mut p = r; for i in r+1..n { if t.setmax(abs(&a[i][c])) { p = i; } } if t == zero { continue } a.swap(r, p); pivots[c] = p; let inv = &one / &a[r][c]; for arj in &mut a[r][c+1..] { *arj = &*arj * &inv } let (ar, a) = a[r..].split_first_mut().unwrap(); for ai in a { let d = ai[c].clone(); for (aij, arj) in ai[c+1..].iter_mut().zip(ar[c+1..].iter()) { *aij -= &d * arj; } } r += 1; if r == n { break } } Decomposed(a, pivots) } } #[derive(Clone, Debug)] pub struct Decomposed(pub Mat, pub Vec); impl Decomposed where T: Clone + Default + PartialOrd + From + SubAssign, for<'b> &'b T: Neg + Mul + Div { pub fn rank(&self) -> usize { let n = self.0.len(); self.1.iter().filter(|&&x| x < n).count() } pub fn det(&self) -> T { let Decomposed(ref a, ref pivots) = *self; assert_eq!(a.len(), a[0].len()); let mut det = T::from(1); for i in 0..a.len() { if pivots[i] == a.len() { return T::default() } det = &det * &a[i][i]; } det } /// Solve Ax=b. x={x[0]+linearspace(x[1..])}. /// O(nr+(1+m-r)r^2 ), r: rank. pub fn solve(&self, mut b: Vec) -> Vec> where T: ::std::fmt::Debug { let Decomposed(ref a, ref pivots) = *self; let (n, m) = (a.len(), a[0].len()); assert_eq!(n, b.len()); let (zero, one) = (T::default(), T::from(1)); let mut id = vec![m; n + 1]; let mut r = 0; for c in 0..m { if pivots[c] == n { continue } b.swap(r, pivots[c]); id[r] = c; r += 1; if r == n { break } } for r in 0..r { let c = id[r]; b[r] = &b[r] / &a[r][c]; for i in r+1..n { let tmp = &a[i][c] * &b[r]; b[i] -= tmp; } } for r in (1..r).rev() { let c = id[r]; for i in 0..r { let tmp = &a[i][c] * &b[r]; b[i] -= tmp; } } if b[r..].iter().any(|v| v != &zero) { return vec![] } let mut x = mat![T::default(); 1 + m - r; m]; let mut k = 0; for j in 0..m { if id[k] == j { x[0][j] = b[k].clone(); k += 1; } else { let mut c: Vec = a[0..k].iter().map(|ai| ai[j].neg()).collect(); for r in (1..k).rev() { for i in 0..r { let tmp = &a[i][id[r]] * &c[r]; c[i] -= tmp; } } for (i, v) in c.into_iter().enumerate() { x[1 + j - k][id[i]] = v; } x[1 + j - k][j] = one.clone(); } } x } } impl Deref for Mat { type Target = Vec>; #[inline] fn deref(&self) -> &Vec> { &self.0 } } impl DerefMut for Mat { #[inline] fn deref_mut(&mut self) -> &mut Vec> { &mut self.0 } } impl<'a, T> Add for &'a Mat where for<'b> &'b T: Add { type Output = Mat; fn add(self, a: &Mat) -> Mat { assert_eq!(self.len(), a.len()); assert_eq!(self[0].len(), a[0].len()); Mat(self.iter().zip(a.iter()).map(|(xs, ys)| xs.iter().zip(ys.iter()).map(|(x, y)| x + y).collect()).collect()) } } impl<'a, T> Sub for &'a Mat where for<'b> &'b T: Sub { type Output = Mat; fn sub(self, a: &Mat) -> Mat { assert_eq!(self.len(), a.len()); assert_eq!(self[0].len(), a[0].len()); Mat(self.iter().zip(a.iter()).map(|(xs, ys)| xs.iter().zip(ys.iter()).map(|(x, y)| x - y).collect()).collect()) } } impl<'a, T> Mul for &'a Mat where T: Clone + Default + AddAssign + PartialEq, for<'b> &'b T: Mul { type Output = Mat; /// 疎行列に対応.n×k行列Aとk×m行列Bの積を O(#nonzero(A) m)で行う. fn mul(self, a: &Mat) -> Mat { assert_eq!(self[0].len(), a.len()); let zero = T::default(); let mut b = mat![zero.clone(); self.len(); a[0].len()]; for (i, si) in self.iter().enumerate() { for (k, sik) in si.iter().enumerate() { if *sik != zero { for (bij, akj) in b[i].iter_mut().zip(a[k].iter()) { *bij += sik * akj; } } } } Mat(b) } } pub struct XorShift { pub x: [u32; 4] } impl XorShift { pub fn new(mut seed: u32) -> XorShift { let mut x = [0; 4]; for i in 0..4 { seed = 1812433253u32.wrapping_mul(seed ^ (seed >> 30)).wrapping_add(i as u32); x[i] = seed; } XorShift { x: x } } pub fn next32(&mut self) -> usize { let t = self.x[0] ^ (self.x[0] << 11); for i in 0..3 { self.x[i] = self.x[i + 1] } self.x[3] = self.x[3] ^ (self.x[3] >> 19) ^ t ^ (t >> 8); self.x[3] as usize } /// [0, n) pub fn next(&mut self, n: usize) -> usize { loop { let t = self.next32(); let r = t % n; if (t - r).checked_add(n).is_some() { return r; } } } } fn readln() -> String { let mut line = String::new(); ::std::io::stdin().read_line(&mut line).unwrap_or_else(|e| panic!("{}", e)); line } #[macro_export] macro_rules! read { ($($t:tt),*; $n:expr) => {{ let stdin = ::std::io::stdin(); let ret = ::std::io::BufRead::lines(stdin.lock()).take($n).map(|line| { let line = line.unwrap(); let mut it = line.split_whitespace(); _read!(it; $($t),*) }).collect::>(); ret }}; ($($t:tt),*) => {{ let line = readln(); let mut it = line.split_whitespace(); _read!(it; $($t),*) }}; } #[macro_export] macro_rules! _read { ($it:ident; [char]) => { _read!($it; String).chars().collect::>() }; ($it:ident; [u8]) => { Vec::from(_read!($it; String).into_bytes()) }; ($it:ident; [$t:ty]) => { $it.map(|s| s.parse::<$t>().unwrap_or_else(|e| panic!("{}", e))).collect::>() }; ($it:ident; $t:ty) => { $it.next().unwrap_or_else(|| panic!("input mismatch")).parse::<$t>().unwrap_or_else(|e| panic!("{}", e)) }; ($it:ident; $($t:tt),+) => { ($(_read!($it; $t)),*) }; }