結果

問題 No.5007 Steiner Space Travel
ユーザー wata_orzwata_orz
提出日時 2022-07-30 17:58:30
言語 Rust
(1.77.0 + proconio)
結果
AC  
実行時間 952 ms / 1,000 ms
コード長 19,411 bytes
コンパイル時間 3,010 ms
実行使用メモリ 3,040 KB
スコア 8,980,319
最終ジャッジ日時 2022-07-30 17:59:25
合計ジャッジ時間 34,043 ms
ジャッジサーバーID
（参考情報） 		judge13 / judge12
純コード判定しない問題か言語
このコードへのチャレンジ
（要ログイン）

テストケース

テストケース表示
入力 結果 実行時間
実行使用メモリ
testcase_00 AC 951 ms
2,996 KB
testcase_01 AC 951 ms
2,976 KB
testcase_02 AC 951 ms
2,948 KB
testcase_03 AC 952 ms
3,040 KB
testcase_04 AC 951 ms
2,948 KB
testcase_05 AC 951 ms
2,980 KB
testcase_06 AC 951 ms
2,964 KB
testcase_07 AC 952 ms
3,012 KB
testcase_08 AC 951 ms
2,968 KB
testcase_09 AC 951 ms
3,004 KB
testcase_10 AC 951 ms
2,968 KB
testcase_11 AC 950 ms
2,988 KB
testcase_12 AC 951 ms
2,972 KB
testcase_13 AC 952 ms
2,980 KB
testcase_14 AC 952 ms
2,956 KB
testcase_15 AC 952 ms
2,992 KB
testcase_16 AC 951 ms
3,036 KB
testcase_17 AC 951 ms
2,984 KB
testcase_18 AC 951 ms
3,012 KB
testcase_19 AC 952 ms
2,960 KB
testcase_20 AC 951 ms
2,860 KB
testcase_21 AC 951 ms
2,884 KB
testcase_22 AC 951 ms
3,008 KB
testcase_23 AC 951 ms
3,016 KB
testcase_24 AC 951 ms
2,876 KB
testcase_25 AC 951 ms
2,988 KB
testcase_26 AC 952 ms
2,976 KB
testcase_27 AC 951 ms
2,976 KB
testcase_28 AC 952 ms
2,960 KB
testcase_29 AC 951 ms
2,992 KB
権限があれば一括ダウンロードができます
コンパイルメッセージ
warning: unnecessary parentheses around block return value
   --> Main.rs:275:4
    |
275 |             (ms - STIME)
    |             ^          ^
    |
    = note: `#[warn(unused_parens)]` on by default
help: remove these parentheses
    |
275 -             (ms - STIME)
275 +             ms - STIME
    | 

warning: unused variable: `rng`
   --> Main.rs:354:63
    |
354 |     pub fn solve(g: &Vec<Vec<i64>>, qs: &Vec<usize>, until: f64, rng: &mut XorShift) -> Vec<usize> {
    |                                                                  ^^^ help: if this is intentional, prefix it with an underscore: `_rng`
    |
    = note: `#[warn(unused_variables)]` on by default

warning: 2 warnings emitted

ソースコード

diff # 

#![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);
			let ty = rng.next(4);
			if out.route.len() > 0 && ty == 0 {
				let mut tmp = vec![];
				for i in 1..out.route.len() {
					if out.route[i - 1] < input.N && out.route[i] < input.N {
						tmp.push((dist2(input.ps[out.route[i - 1]], input.ps[out.route[i]]), i));
					}
				}
				if tmp.len() == 0 {
					continue;
				}
				tmp.sort();
				tmp.reverse();
				let i = tmp[rng.next(tmp.len().min(5))].1;
				ps[input.N + p].0 = (input.ps[out.route[i - 1]].0 + input.ps[out.route[i]].0) / 2;
				ps[input.N + p].1 = (input.ps[out.route[i - 1]].1 + input.ps[out.route[i]].1) / 2;
			} else if ty < 3 {
				ps[input.N + p].0 += rng.next(201) as i32 - 100;
				ps[input.N + p].1 += rng.next(201) as i32 - 100;
			} else {
				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<usize>,
}

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<T> 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<Vec<i64>>, ps: &Vec<usize>) -> 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<i64>>) -> Vec<usize> {
		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<usize>, idx: &mut Vec<usize>, 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<Vec<i64>>, qs: &Vec<usize>, until: f64, rng: &mut XorShift) -> Vec<usize> {
		// 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<T>(pub Vec<Vec<T>>);

impl<T> Mat<T> where T: Clone + Default {
	pub fn row(a: Vec<T>) -> Mat<T> { Mat(vec![a]) }
	pub fn col(a: Vec<T>) -> Mat<T> { Mat(a.into_iter().map(|v| vec![v]).collect()) }
	pub fn to_vec(&self) -> Vec<T> {
		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<T> {
		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<T> where T: AddAssign + PartialEq + From<u8>, for<'b> &'b T: Mul<Output = T> {
		let n = self.len();
		assert_eq!(n, self[0].len());
		let mut res: Mat<T> = 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<T> where T: Clone + Default + PartialOrd + From<u8> + SubAssign, for<'b> &'b T: Neg<Output = T> + Mul<Output = T> + Div<Output = T> {
		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<T>(pub Mat<T>, pub Vec<usize>);

impl<T> Decomposed<T> where T: Clone + Default + PartialOrd + From<u8> + SubAssign, for<'b> &'b T: Neg<Output = T> + Mul<Output = T> + Div<Output = T> {
	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<T>) -> Vec<Vec<T>> 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<T> = 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<T> Deref for Mat<T> {
	type Target = Vec<Vec<T>>;
	#[inline]
	fn deref(&self) -> &Vec<Vec<T>> { &self.0 }
}

impl<T> DerefMut for Mat<T> {
	#[inline]
	fn deref_mut(&mut self) -> &mut Vec<Vec<T>> { &mut self.0 }
}

impl<'a, T> Add for &'a Mat<T> where for<'b> &'b T: Add<Output = T> {
	type Output = Mat<T>;
	fn add(self, a: &Mat<T>) -> Mat<T> {
		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<T> where for<'b> &'b T: Sub<Output = T> {
	type Output = Mat<T>;
	fn sub(self, a: &Mat<T>) -> Mat<T> {
		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<T> where T: Clone + Default + AddAssign + PartialEq, for<'b> &'b T: Mul<Output = T> {
	type Output = Mat<T>;
	/// 疎行列に対応．n×k行列Aとk×m行列Bの積を O(#nonzero(A) m)で行う．
	fn mul(self, a: &Mat<T>) -> Mat<T> {
		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::<Vec<_>>();
		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::<Vec<_>>()
	};
	($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::<Vec<_>>()
	};
	($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)),*)
	};
}
0