// codesnip-guard: BlackBoxMatrix
#[cfg_attr(any(),rustfmt::skip)]pub use self::black_box_matrix::{BlackBoxMatrix,BlackBoxMatrixImpl,SparseMatrix};#[cfg_attr(any(),rustfmt::skip)]mod black_box_matrix{use super::{Field,Invertible,Matrix,SemiRing};use std::{fmt::{self,Debug},marker::PhantomData};pub trait BlackBoxMatrix<R>where R:SemiRing{fn apply(&self,v:&[R::T])->Vec<R::T>;fn shape(&self)->(usize,usize);}impl<R>BlackBoxMatrix<R>for Matrix<R>where R:SemiRing{fn apply(&self,v:&[R::T])->Vec<R::T>{assert_eq!(self.shape.1,v.len());self.data.iter().map(|row|R::dot_product(row,v)).collect()}fn shape(&self)->(usize,usize){self.shape}}pub struct SparseMatrix<R>where R:SemiRing{shape:(usize,usize),nonzero:Vec<(usize,usize,R::T)>}impl<R>Debug for SparseMatrix<R>where R:SemiRing<T:Debug>{fn fmt(&self,f:&mut fmt::Formatter<'_>)->fmt::Result{f.debug_struct("SparseMatrix").field("shape",&self.shape).field("nonzero",&self.nonzero).finish()}}impl<R>Clone for SparseMatrix<R>where R:SemiRing{fn clone(&self)->Self{Self{shape:self.shape,nonzero:self.nonzero.clone()}}}impl<R>SparseMatrix<R>where R:SemiRing{pub fn new(shape:(usize,usize))->Self{Self{shape,nonzero:vec![]}}pub fn new_with<F>(shape:(usize,usize),f:F)->Self where R:SemiRing<T:PartialEq>,F:Fn(usize,usize)->R::T{let mut nonzero=vec![];for i in 0..shape.0{for j in 0..shape.1{let v=f(i,j);if!R::is_zero(&v){nonzero.push((i,j,v));}}}Self{shape,nonzero}}pub fn from_nonzero(shape:(usize,usize),nonzero:Vec<(usize,usize,R::T)>)->Self{Self{shape,nonzero}}}impl<R>SparseMatrix<R>where R:Field<T:PartialEq,Additive:Invertible,Multiplicative:Invertible>{pub fn determinant(&self)->R::T{assert_eq!(self.shape.0,self.shape.1);let n=self.shape.0;let mut columns=vec![Vec::<(usize,R::T)>::new();n];for&(i,j,ref value)in&self.nonzero{columns[j].push((i,value.clone()));}let mut degrees=vec![0;n];for column in&mut columns{column.sort_unstable_by_key(|&(i,_)|i);let mut merged:Vec<(usize,R::T)>=Vec::with_capacity(column.len());for(i,value)in column.drain(..){if let Some((_,x))=merged.last_mut().filter(|(last,_)|*last==i){R::add_assign(x,&value);}else{merged.push((i,value));}}merged.retain(|(i,value)|{if R::is_zero(value){false}else{degrees[*i]+=1;true}});*column=merged;}let mut order:Vec<_>=(0..n).collect();order.sort_unstable_by_key(|&j|columns[j].len());let mut lower:Vec<Vec<(usize,R::T)>>=Vec::with_capacity(n);let mut pivots:Vec<Option<usize>>=vec![None;n];let mut x=vec![R::zero();n];let mut seen=vec![0;n];let mut stack=Vec::new();let mut support=Vec::new();let mut determinant=R::one();for(k,&j)in order.iter().enumerate(){support.clear();for&(i,_)in&columns[j]{if seen[i]==k+1{continue;}seen[i]=k+1;x[i]=R::zero();stack.push((i,0));while let Some((i,next))=stack.last_mut(){if let Some(pivot)=pivots[*i]&&*next<lower[pivot].len(){let row=lower[pivot][*next].0;*next+=1;if seen[row]!=k+1{seen[row]=k+1;x[row]=R::zero();stack.push((row,0));}continue;}support.push(*i);stack.pop();}}for&(i,ref value)in&columns[j]{x[i]=value.clone();}let mut pivot=None;for&i in support.iter().rev(){if let Some(p)=pivots[i]{let factor=R::neg(&x[i]);for&(row,ref value)in&lower[p]{R::add_assign(&mut x[row],&R::mul(&factor,value));}}else if!R::is_zero(&x[i])&&pivot.is_none_or(|p|degrees[i]<degrees[p]){pivot=Some(i);}}let Some(pivot)=pivot else{return R::zero()};R::mul_assign(&mut determinant,&x[pivot]);let inv=R::inv(&x[pivot]);pivots[pivot]=Some(k);lower.push(support.iter().filter(|&&i|pivots[i].is_none()&&!R::is_zero(&x[i])).map(|&i|(i,R::mul(&x[i],&inv))).collect());}for mut permutation in[order,pivots.into_iter().map(Option::unwrap).collect()]{for i in 0..n{while permutation[i]!=i{let j=permutation[i];permutation.swap(i,j);determinant=R::neg(&determinant);}}}determinant}}impl<R>From<Matrix<R>>for SparseMatrix<R>where R:SemiRing<T:PartialEq>{fn from(mat:Matrix<R>)->Self{let mut nonzero=vec![];for i in 0..mat.shape.0{for j in 0..mat.shape.1{let v=mat[(i,j)].clone();if!R::is_zero(&v){nonzero.push((i,j,v));}}}Self{shape:mat.shape,nonzero}}}impl<R>From<SparseMatrix<R>>for Matrix<R>where R:SemiRing{fn from(smat:SparseMatrix<R>)->Self{let mut mat=Matrix::zeros(smat.shape);for&(i,j,ref v)in&smat.nonzero{R::add_assign(&mut mat[(i,j)],v);}mat}}impl<R>BlackBoxMatrix<R>for SparseMatrix<R>where R:SemiRing{fn apply(&self,v:&[R::T])->Vec<R::T>{assert_eq!(self.shape.1,v.len());let mut res=vec![R::zero();self.shape.0];for&(i,j,ref val)in&self.nonzero{R::add_assign(&mut res[i],&R::mul(val,&v[j]));}res}fn shape(&self)->(usize,usize){self.shape}}pub struct BlackBoxMatrixImpl<R,F>{shape:(usize,usize),apply_fn:F,_marker:PhantomData<fn()->R>}impl<R,F>Debug for BlackBoxMatrixImpl<R,F>where F:Debug{fn fmt(&self,f:&mut fmt::Formatter<'_>)->fmt::Result{f.debug_struct("BlackBoxMatrixImpl").field("shape",&self.shape).field("apply_fn",&self.apply_fn).finish()}}impl<R,F>Clone for BlackBoxMatrixImpl<R,F>where F:Clone{fn clone(&self)->Self{Self{shape:self.shape,apply_fn:self.apply_fn.clone(),_marker:PhantomData}}}impl<R,F>BlackBoxMatrixImpl<R,F>{pub fn new(shape:(usize,usize),apply_fn:F)->Self{Self{shape,apply_fn,_marker:PhantomData}}}impl<R,F>BlackBoxMatrix<R>for BlackBoxMatrixImpl<R,F>where R:SemiRing,F:Fn(&[R::T])->Vec<R::T>{fn apply(&self,v:&[R::T])->Vec<R::T>{assert_eq!(self.shape.1,v.len());(self.apply_fn)(v)}fn shape(&self)->(usize,usize){self.shape}}}
// codesnip-guard: AdditiveOperation
#[cfg_attr(any(),rustfmt::skip)]pub use self::additive_operation_impl::AdditiveOperation;#[cfg_attr(any(),rustfmt::skip)]mod additive_operation_impl{use super::*;use std::{marker::PhantomData,ops::{Add,Neg,Sub}};#[doc=" $+$"]pub struct AdditiveOperation<T>where T:Clone+Zero+Add<Output=T>{_marker:PhantomData<fn()->T>}impl<T>Magma for AdditiveOperation<T>where T:Clone+Zero+Add<Output=T>{type T=T;#[inline]fn operate(x:&Self::T,y:&Self::T)->Self::T{x.clone()+y.clone()}}impl<T>Unital for AdditiveOperation<T>where T:Clone+Zero+Add<Output=T>{#[inline]fn unit()->Self::T{Zero::zero()}}impl<T>Associative for AdditiveOperation<T>where T:Clone+Zero+Add<Output=T>{}impl<T>Commutative for AdditiveOperation<T>where T:Clone+Zero+Add<Output=T>{}impl<T>Invertible for AdditiveOperation<T>where T:Clone+Zero+Add<Output=T>+Sub<Output=T>+Neg<Output=T>{#[inline]fn inverse(x:&Self::T)->Self::T{-x.clone()}#[inline]fn rinv_operate(x:&Self::T,y:&Self::T)->Self::T{x.clone()-y.clone()}}}
// codesnip-guard: Matrix
#[cfg_attr(any(),rustfmt::skip)]pub use self::matrix::Matrix;#[cfg_attr(any(),rustfmt::skip)]mod matrix{use super::{Field,Invertible,Ring,SemiRing};use std::{fmt::{self,Debug},marker::PhantomData,ops::{Add,AddAssign,Index,IndexMut,Mul,MulAssign,Neg,Sub,SubAssign}};pub struct Matrix<R>where R:SemiRing{pub shape:(usize,usize),pub data:Vec<Vec<R::T>>,_marker:PhantomData<fn()->R>}impl<R>Debug for Matrix<R>where R:SemiRing<T:Debug>{fn fmt(&self,f:&mut fmt::Formatter<'_>)->fmt::Result{f.debug_struct("Matrix").field("shape",&self.shape).field("data",&self.data).field("_marker",&self._marker).finish()}}impl<R>Clone for Matrix<R>where R:SemiRing{fn clone(&self)->Self{Self{shape:self.shape,data:self.data.clone(),_marker:self._marker}}}impl<R>PartialEq for Matrix<R>where R:SemiRing<T:PartialEq>{fn eq(&self,other:&Self)->bool{self.shape==other.shape&&self.data==other.data}}impl<R>Eq for Matrix<R>where R:SemiRing<T:Eq>{}impl<R>Matrix<R>where R:SemiRing{pub fn new(shape:(usize,usize),z:R::T)->Self{Self{shape,data:vec![vec![z;shape.1];shape.0],_marker:PhantomData}}pub fn from_vec(data:Vec<Vec<R::T>>)->Self{let shape=(data.len(),data.first().map(Vec::len).unwrap_or_default());assert!(data.iter().all(|r|r.len()==shape.1));Self{shape,data,_marker:PhantomData}}pub fn new_with(shape:(usize,usize),mut f:impl FnMut(usize,usize)->R::T)->Self{let data=(0..shape.0).map(|i|(0..shape.1).map(|j|f(i,j)).collect()).collect();Self{shape,data,_marker:PhantomData}}pub fn zeros(shape:(usize,usize))->Self{Self{shape,data:vec![vec![R::zero();shape.1];shape.0],_marker:PhantomData}}pub fn eye(shape:(usize,usize))->Self{let mut data=vec![vec![R::zero();shape.1];shape.0];for(i,d)in data.iter_mut().enumerate().take(shape.1){d[i]=R::one();}Self{shape,data,_marker:PhantomData}}pub fn transpose(&self)->Self{Self::new_with((self.shape.1,self.shape.0),|i,j|self[j][i].clone())}pub fn map<S,F>(&self,mut f:F)->Matrix<S>where S:SemiRing,F:FnMut(&R::T)->S::T{Matrix::<S>::new_with(self.shape,|i,j|f(&self[i][j]))}pub fn add_row_with(&mut self,mut f:impl FnMut(usize,usize)->R::T){self.data.push((0..self.shape.1).map(|j|f(self.shape.0,j)).collect());self.shape.0+=1;}pub fn add_col_with(&mut self,mut f:impl FnMut(usize,usize)->R::T){for i in 0..self.shape.0{self.data[i].push(f(i,self.shape.1));}self.shape.1+=1;}pub fn pairwise_assign<F>(&mut self,other:&Self,mut f:F)where F:FnMut(&mut R::T,&R::T){assert_eq!(self.shape,other.shape);for i in 0..self.shape.0{for j in 0..self.shape.1{f(&mut self[i][j],&other[i][j]);}}}}#[derive(Debug)]pub struct SystemOfLinearEquationsSolution<R>where R:Field<Additive:Invertible,Multiplicative:Invertible>{pub particular:Vec<R::T>,pub basis:Vec<Vec<R::T>>}impl<R>Matrix<R>where R:Field<T:PartialEq,Additive:Invertible,Multiplicative:Invertible>{fn eliminate<const DETERMINANT:bool>(&mut self)->(usize,R::T){let(n,m)=self.shape;let mut rank=0;let mut determinant=R::one();let mut negative=false;for first in(0..m).step_by(64){if rank==n{break;}let end=(first+64).min(m);let start=rank;let mut pivots=Vec::new();let mut panel=vec![vec![R::zero();end-first];end-first];for col in first..end{if panel[col-first][..col-first].iter().any(|x|!R::is_zero(x)){for row in&mut self.data[rank..]{let value=R::dot_product(&row[first..col],&panel[col-first][..col-first]);R::sub_assign(&mut row[col],&value);}}let Some(pivot)=(rank..n).find(|&i|!R::is_zero(&self[i][col]))else{continue;};if pivot!=rank{self.data.swap(rank,pivot);negative=!negative;}if DETERMINANT{R::mul_assign(&mut determinant,&self[rank][col]);}let inv=R::inv(&self[rank][col]);let row=&mut self.data[rank];for c in col+1..end{let value=R::dot_product(&row[first..col],&panel[c-first][..col-first]);R::sub_assign(&mut row[c],&value);panel[c-first][col-first]=row[c].clone();}for row in&mut self.data[rank+1..]{R::mul_assign(&mut row[col],&inv);}pivots.push(col);rank+=1;if rank==n{break;}}for i in start..if rank-start<32{n}else{rank}{let(upper,lower)=self.data.split_at_mut(i);let row=&mut lower[0];for(j,&col)in pivots[..(i-start).min(pivots.len())].iter().enumerate(){if R::is_zero(&row[col]){continue;}let factor=R::neg(&row[col]);R::add_scaled_assign(&mut row[end..],&upper[start+j][end..],&factor);}}if rank<n&&end<m&&rank-start>=32&&self.data[rank..].iter().any(|row|pivots.iter().any(|&col|!R::is_zero(&row[col]))){let lower=Self::new_with((n-rank,rank-start),|i,j|{R::neg(&self[rank+i][pivots[j]])});let upper=Self::new_with((rank-start,m-end),|i,j|{self[start+i][end+j].clone()});let update=&lower*&upper;for(row,update)in self.data[rank..].iter_mut().zip(update.data){for(x,y)in row[end..].iter_mut().zip(update){R::add_assign(x,&y);}}}for(i,&col)in pivots.iter().enumerate(){for row in&mut self.data[start+i+1..]{row[col]=R::zero();}}if DETERMINANT&&rank<end{return(rank,R::zero());}}if DETERMINANT&&negative{determinant=R::neg(&determinant);}(rank,determinant)}#[doc=" f: (row, pivot_row, col)"]pub fn row_reduction_with<F>(&mut self,normalize:bool,mut f:F)where F:FnMut(usize,usize,usize){let(n,m)=self.shape;let mut c=0;for r in 0..n{loop{if c>=m{return;}if let Some(pivot)=(r..n).find(|&p|!R::is_zero(&self[p][c])){f(r,pivot,c);self.data.swap(r,pivot);break;};c+=1;}let d=R::inv(&self[r][c]);if normalize{for value in&mut self[r][c..m]{R::mul_assign(value,&d);}}for i in(0..n).filter(|&i|i!=r){let mut e=self[i][c].clone();if!normalize{R::mul_assign(&mut e,&d);}for j in c..m{let e=R::mul(&e,&self[r][j]);R::sub_assign(&mut self[i][j],&e);}}c+=1;}}pub fn row_reduction(&mut self,normalize:bool){self.row_reduction_with(normalize,|_,_,_|{});}pub fn rank(&mut self)->usize{self.eliminate::<false>().0}pub fn determinant(&mut self)->R::T{assert_eq!(self.shape.0,self.shape.1);self.eliminate::<true>().1}pub fn solve_system_of_linear_equations(&self,b:&[R::T])->Option<SystemOfLinearEquationsSolution<R>>{assert_eq!(self.shape.0,b.len());let m=self.shape.1;let mut a=Self::new_with((self.shape.0,m+1),|i,j|{if j==m{b[i].clone()}else{self[i][j].clone()}});let rank=a.eliminate::<false>().0;let mut pivots=Vec::with_capacity(rank);let mut b=Vec::with_capacity(rank);for row in&a.data[..rank]{let c=row.iter().position(|x|!R::is_zero(x)).unwrap();if c==m{return None;}pivots.push(c);b.push(row[m].clone());}let mut free=Vec::with_capacity(m-rank);let mut pivot=0;for c in 0..m{if pivot<rank&&pivots[pivot]==c{pivot+=1;}else{free.push(c);}}let mut coefficients:Vec<Vec<_>>=(0..rank).map(|i|free.iter().map(|&c|a[i][c].clone()).collect()).collect();for k in(0..rank).rev(){let c=pivots[k];let inv=R::inv(&a[k][c]);R::mul_assign(&mut b[k],&inv);let pivot_b=b[k].clone();let(upper,lower)=coefficients.split_at_mut(k);let pivot_coefficients=&mut lower[0];for x in pivot_coefficients.iter_mut(){R::mul_assign(x,&inv);}for((row,value),coefficients)in a.data[..k].iter_mut().zip(&mut b[..k]).zip(upper){if R::is_zero(&row[c]){continue;}let factor=row[c].clone();row[c]=R::zero();R::sub_assign(value,&R::mul(&factor,&pivot_b));R::add_scaled_assign(coefficients,pivot_coefficients,&R::neg(&factor));}}let mut particular=vec![R::zero();m];for i in 0..rank{particular[pivots[i]]=b[i].clone();}let mut basis=Vec::with_capacity(free.len());for(j,&c)in free.iter().enumerate(){let mut vector=vec![R::zero();m];vector[c]=R::one();for i in 0..rank{vector[pivots[i]]=R::neg(&coefficients[i][j]);}basis.push(vector);}Some(SystemOfLinearEquationsSolution{particular,basis})}pub fn inverse(&self)->Option<Matrix<R>>{assert_eq!(self.shape.0,self.shape.1);let n=self.shape.0;if n>=64{let m=n/2;let a=Self::new_with((m,m),|i,j|self[i][j].clone());if let Some(mut ai)=a.inverse(){let b=Self::new_with((m,n-m),|i,j|self[i][j+m].clone());let c=Self::new_with((n-m,m),|i,j|self[i+m][j].clone());let mut d=Self::new_with((n-m,n-m),|i,j|self[i+m][j+m].clone());let u=&ai*&b;let v=&c*&ai;d-=&v*&b;let di=d.inverse()?;let r=&u*&di;let t=&di*&v;ai+=&r*&v;let mut inverse=Self::zeros((n,n));for i in 0..m{inverse[i][..m].clone_from_slice(&ai[i]);for(x,y)in inverse[i][m..].iter_mut().zip(&r[i]){*x=R::neg(y);}}for i in m..n{for(x,y)in inverse[i][..m].iter_mut().zip(&t[i-m]){*x=R::neg(y);}inverse[i][m..].clone_from_slice(&di[i-m]);}return Some(inverse);}}let mut a=self.clone();let mut inverse=Self::eye((n,n));let mut ranges:Vec<_>=(0..n).map(|i|(i,i+1)).collect();for r in 0..n{let pivot=(r..n).find(|&i|!R::is_zero(&a[i][r]))?;a.data.swap(r,pivot);inverse.data.swap(r,pivot);ranges.swap(r,pivot);let d=R::inv(&a[r][r]);for x in&mut a[r][r..]{R::mul_assign(x,&d);}let(left,right)=ranges[r];for x in&mut inverse[r][left..right]{R::mul_assign(x,&d);}let(a_upper,a_lower)=a.data.split_at_mut(r+1);let pivot_a=&a_upper[r];let(inverse_upper,inverse_lower)=inverse.data.split_at_mut(r+1);let pivot_inverse=&inverse_upper[r];let(ranges_upper,ranges_lower)=ranges.split_at_mut(r+1);let(left,right)=ranges_upper[r];for((a,inverse),range)in a_lower.iter_mut().zip(inverse_lower).zip(ranges_lower){if R::is_zero(&a[r]){continue;}let e=a[r].clone();a[r]=R::zero();R::add_scaled_assign(&mut a[(r+1)..],&pivot_a[(r+1)..],&R::neg(&e));R::add_scaled_assign(&mut inverse[left..right],&pivot_inverse[left..right],&R::neg(&e));range.0=range.0 .min(left);range.1=range.1 .max(right);}}for r in(0..n).rev(){let(left,right)=ranges[r];let(inverse_upper,inverse_lower)=inverse.data.split_at_mut(r);let pivot_inverse=&inverse_lower[0];let(ranges_upper,_)=ranges.split_at_mut(r);for((a,inverse),range)in a.data[..r].iter_mut().zip(inverse_upper).zip(ranges_upper){if R::is_zero(&a[r]){continue;}let e=a[r].clone();a[r]=R::zero();R::add_scaled_assign(&mut inverse[left..right],&pivot_inverse[left..right],&R::neg(&e));range.0=range.0 .min(left);range.1=range.1 .max(right);}}Some(inverse)}pub fn characteristic_polynomial(&mut self)->Vec<R::T>{let n=self.shape.0;if n==0{return vec![R::one()];}assert!(self.data.iter().all(|a|a.len()==n));for j in 0..(n-1){if let Some(x)=((j+1)..n).find(|&x|!R::is_zero(&self[x][j])){self.data.swap(j+1,x);self.data.iter_mut().for_each(|a|a.swap(j+1,x));let inv=R::inv(&self[j+1][j]);let mut v=vec![];let src=std::mem::take(&mut self[j+1]);for a in self.data[(j+2)..].iter_mut(){let mul=R::mul(&a[j],&inv);R::add_scaled_assign(&mut a[j..],&src[j..],&R::neg(&mul));v.push(mul);}self[j+1]=src;for a in self.data.iter_mut(){let v=R::dot_product(&a[(j+2)..],&v);R::add_assign(&mut a[j+1],&v);}}}let mut dp:Vec<Vec<R::T>>=(0..=n).map(|i|Vec::with_capacity(n+1-i)).collect();dp[0].push(R::one());for i in 0..n{let mut c=vec![R::zero();i+1];c[i]=R::neg(&self[i][i]);let mut mul=R::one();for j in(0..i).rev(){mul=R::mul(&mul,&self[j+1][j]);c[j]=R::neg(&R::mul(&mul,&self[j][i]));}for k in(0..=i).rev(){let mut value=R::dot_product(&dp[k],&c[k..]);if k>0{R::add_assign(&mut value,dp[k-1].last().unwrap());}dp[k].push(value);}dp[i+1].push(R::one());}dp.into_iter().map(|mut c|c.pop().unwrap()).collect()}}impl<R>Index<usize>for Matrix<R>where R:SemiRing{type Output=Vec<R::T>;fn index(&self,index:usize)->&Self::Output{&self.data[index]}}impl<R>IndexMut<usize>for Matrix<R>where R:SemiRing{fn index_mut(&mut self,index:usize)->&mut Self::Output{&mut self.data[index]}}impl<R>Index<(usize,usize)>for Matrix<R>where R:SemiRing{type Output=R::T;fn index(&self,index:(usize,usize))->&Self::Output{&self.data[index.0][index.1]}}impl<R>IndexMut<(usize,usize)>for Matrix<R>where R:SemiRing{fn index_mut(&mut self,index:(usize,usize))->&mut Self::Output{&mut self.data[index.0][index.1]}}macro_rules!impl_matrix_pairwise_binop{($imp:ident,$method:ident,$imp_assign:ident,$method_assign:ident$(where[$($clauses:tt)*])?)=>{impl<R>$imp_assign for Matrix<R>where R:SemiRing,$($($clauses)*)?{fn$method_assign(&mut self,rhs:Self){self.pairwise_assign(&rhs,|a,b|R::$method_assign(a,b));}}impl<R>$imp_assign<&Matrix<R>>for Matrix<R>where R:SemiRing,$($($clauses)*)?{fn$method_assign(&mut self,rhs:&Self){self.pairwise_assign(rhs,|a,b|R::$method_assign(a,b));}}impl<R>$imp for Matrix<R>where R:SemiRing,$($($clauses)*)?{type Output=Matrix<R>;fn$method(mut self,rhs:Self)->Self::Output{self.$method_assign(rhs);self}}impl<R>$imp<&Matrix<R>>for Matrix<R>where R:SemiRing,$($($clauses)*)?{type Output=Matrix<R>;fn$method(mut self,rhs:&Self)->Self::Output{self.$method_assign(rhs);self}}impl<R>$imp<Matrix<R>>for&Matrix<R>where R:SemiRing,$($($clauses)*)?{type Output=Matrix<R>;fn$method(self,mut rhs:Matrix<R>)->Self::Output{rhs.pairwise_assign(self,|a,b|*a=R::$method(b,a));rhs}}impl<R>$imp<&Matrix<R>>for&Matrix<R>where R:SemiRing,$($($clauses)*)?{type Output=Matrix<R>;fn$method(self,rhs:&Matrix<R>)->Self::Output{let mut this=self.clone();this.$method_assign(rhs);this}}};}impl_matrix_pairwise_binop!(Add,add,AddAssign,add_assign);impl_matrix_pairwise_binop!(Sub,sub,SubAssign,sub_assign where[R:SemiRing<Additive:Invertible>]);impl<R>Mul for Matrix<R>where R:SemiRing{type Output=Matrix<R>;fn mul(self,rhs:Self)->Self::Output{(&self).mul(&rhs)}}impl<R>Mul<&Matrix<R>>for Matrix<R>where R:SemiRing{type Output=Matrix<R>;fn mul(self,rhs:&Matrix<R>)->Self::Output{(&self).mul(rhs)}}impl<R>Mul<Matrix<R>>for&Matrix<R>where R:SemiRing{type Output=Matrix<R>;fn mul(self,rhs:Matrix<R>)->Self::Output{self.mul(&rhs)}}impl<R>Mul<&Matrix<R>>for&Matrix<R>where R:SemiRing{type Output=Matrix<R>;fn mul(self,rhs:&Matrix<R>)->Self::Output{assert_eq!(self.shape.1,rhs.shape.0);if let Some(data)=R::try_matrix_product(&self.data,&rhs.data){return Matrix::from_vec(data);}let rhs=rhs.transpose();Matrix::new_with((self.shape.0,rhs.shape.0),|i,j|{R::dot_product(&self[i],&rhs[j])})}}fn strassen_rec<R:Ring>(a:&[R::T],b:&[R::T],c:&mut[R::T],shape:(usize,usize,usize),stride_a:usize,stride_b:usize){let(n,m,p)=shape;fn add_block<R:Ring>(a:&[R::T],b:&[R::T],out:&mut[R::T],n:usize,stride_a:usize,stride_b:usize){for((a,b),c)in a.chunks(stride_a).zip(b.chunks(stride_b)).zip(out.chunks_exact_mut(n)){for((a,b),c)in a.iter().zip(b.iter()).zip(c.iter_mut()){*c=R::add(a,b);}}}fn sub_block<R:Ring>(a:&[R::T],b:&[R::T],out:&mut[R::T],n:usize,stride_a:usize,stride_b:usize){for((a,b),c)in a.chunks(stride_a).zip(b.chunks(stride_b)).zip(out.chunks_exact_mut(n)){for((a,b),c)in a.iter().zip(b.iter()).zip(c.iter_mut()){*c=R::sub(a,b);}}}if n.min(m).min(p)<=128{let transposed:Vec<_>=(0..p).flat_map(|j|(0..m).map(move|i|b[i*stride_b+j].clone())).collect();for(a,c)in a.chunks(stride_a).zip(c.chunks_exact_mut(p)){for(b,c)in transposed.chunks_exact(m).zip(c){*c=R::dot_product(&a[..m],b);}}return;}let(h,k,w)=(n/2,m/2,p/2);let a11=0;let a12=k;let a21=h*stride_a;let a22=a21+k;let b11=0;let b12=w;let b21=k*stride_b;let b22=b21+w;let block=h*w;let mut buf=vec![R::zero();h*k+k*w+block*7];let(s1,rest)=buf.split_at_mut(h*k);let(s2,m_buf)=rest.split_at_mut(k*w);let(m1,rest)=m_buf.split_at_mut(block);let(m2,rest)=rest.split_at_mut(block);let(m3,rest)=rest.split_at_mut(block);let(m4,rest)=rest.split_at_mut(block);let(m5,rest)=rest.split_at_mut(block);let(m6,m7)=rest.split_at_mut(block);add_block::<R>(&a[a11..],&a[a22..],s1,k,stride_a,stride_a);add_block::<R>(&b[b11..],&b[b22..],s2,w,stride_b,stride_b);strassen_rec::<R>(s1,s2,m1,(h,k,w),k,w);add_block::<R>(&a[a21..],&a[a22..],s1,k,stride_a,stride_a);strassen_rec::<R>(s1,&b[b11..],m2,(h,k,w),k,stride_b);sub_block::<R>(&b[b12..],&b[b22..],s2,w,stride_b,stride_b);strassen_rec::<R>(&a[a11..],s2,m3,(h,k,w),stride_a,w);sub_block::<R>(&b[b21..],&b[b11..],s2,w,stride_b,stride_b);strassen_rec::<R>(&a[a22..],s2,m4,(h,k,w),stride_a,w);add_block::<R>(&a[a11..],&a[a12..],s1,k,stride_a,stride_a);strassen_rec::<R>(s1,&b[b22..],m5,(h,k,w),k,stride_b);sub_block::<R>(&a[a21..],&a[a11..],s1,k,stride_a,stride_a);add_block::<R>(&b[b11..],&b[b12..],s2,w,stride_b,stride_b);strassen_rec::<R>(s1,s2,m6,(h,k,w),k,w);sub_block::<R>(&a[a12..],&a[a22..],s1,k,stride_a,stride_a);add_block::<R>(&b[b21..],&b[b22..],s2,w,stride_b,stride_b);strassen_rec::<R>(s1,s2,m7,(h,k,w),k,w);let c11=0;let c12=w;let c21=h*p;let c22=c21+w;for((((m1,m4),m5),m7),c)in m1.iter().zip(m4.iter()).zip(m5.iter()).zip(m7.iter()).zip(c[c11..].chunks_mut(p).flat_map(|c|c.iter_mut().take(w))){*c=R::add(m1,m4);R::sub_assign(c,m5);R::add_assign(c,m7);}for((m3,m5),c)in m3.iter().zip(m5.iter()).zip(c[c12..].chunks_mut(p).flat_map(|c|c.iter_mut().take(w))){*c=R::add(m3,m5);}for((m2,m4),c)in m2.iter().zip(m4.iter()).zip(c[c21..].chunks_mut(p).flat_map(|c|c.iter_mut().take(w))){*c=R::add(m2,m4);}for((((m1,m2),m3),m6),c)in m1.iter().zip(m2.iter()).zip(m3.iter()).zip(m6.iter()).zip(c[c22..].chunks_mut(p).flat_map(|c|c.iter_mut().take(w))){*c=R::sub(m1,m2);R::add_assign(c,m3);R::add_assign(c,m6);}}impl<R>Matrix<R>where R:Ring{pub fn mul_strassen(&self,rhs:&Matrix<R>)->Matrix<R>{assert_eq!(self.shape.1,rhs.shape.0);if let Some(data)=R::try_matrix_product(&self.data,&rhs.data){return Matrix::from_vec(data);}let(n,m)=self.shape;let p=rhs.shape.1;if n==0||m==0||p==0{return Matrix::zeros((n,p));}let split=n.min(m).min(p).div_ceil(128).next_power_of_two();if split<=2{return self*rhs;}let rows=n.div_ceil(split)*split;let inner=m.div_ceil(split)*split;let cols=p.div_ceil(split)*split;let mut a=vec![R::zero();rows*inner];for(a,data)in a.chunks_exact_mut(inner).zip(&self.data){a[..m].clone_from_slice(data);}let mut b=vec![R::zero();inner*cols];for(b,data)in b.chunks_exact_mut(cols).zip(&rhs.data){b[..p].clone_from_slice(data);}let mut c=vec![R::zero();rows*cols];strassen_rec::<R>(&a,&b,&mut c,(rows,inner,cols),inner,cols);let mut res=Matrix::zeros((n,p));for(data,c)in res.data.iter_mut().zip(c.chunks_exact(cols)){data.clone_from_slice(&c[..p]);}res}}impl<R>MulAssign<&R::T>for Matrix<R>where R:SemiRing{fn mul_assign(&mut self,rhs:&R::T){for i in 0..self.shape.0{for j in 0..self.shape.1{R::mul_assign(&mut self[(i,j)],rhs);}}}}impl<R>Neg for Matrix<R>where R:SemiRing<Additive:Invertible>{type Output=Self;fn neg(self)->Self::Output{self.map(|x|R::neg(x))}}impl<R>Neg for&Matrix<R>where R:SemiRing<Additive:Invertible>{type Output=Matrix<R>;fn neg(self)->Self::Output{self.map(|x|R::neg(x))}}impl<R>Matrix<R>where R:SemiRing{pub fn pow(self,mut n:usize)->Self{assert_eq!(self.shape.0,self.shape.1);let mut res=Matrix::eye(self.shape);let mut x=self;while n>0{if n&1==1{res=&res*&x;}n>>=1;if n!=0{x=&x*&x;}}res}}impl<R>Matrix<R>where R:Ring{pub fn pow_strassen(self,mut n:usize)->Self{assert_eq!(self.shape.0,self.shape.1);let mut res=Matrix::eye(self.shape);let mut x=self;while n>0{if n&1==1{res=res.mul_strassen(&x);}n>>=1;if n!=0{x=x.mul_strassen(&x);}}res}}}
// codesnip-guard: MultiplicativeOperation
#[cfg_attr(any(),rustfmt::skip)]pub use self::multiplicative_operation_impl::MultiplicativeOperation;#[cfg_attr(any(),rustfmt::skip)]mod multiplicative_operation_impl{use super::*;use std::{marker::PhantomData,ops::{Div,Mul}};#[doc=" $\\times$"]pub struct MultiplicativeOperation<T>where T:Clone+One+Mul<Output=T>{_marker:PhantomData<fn()->T>}impl<T>Magma for MultiplicativeOperation<T>where T:Clone+One+Mul<Output=T>{type T=T;#[inline]fn operate(x:&Self::T,y:&Self::T)->Self::T{x.clone()*y.clone()}}impl<T>Unital for MultiplicativeOperation<T>where T:Clone+One+Mul<Output=T>{#[inline]fn unit()->Self::T{One::one()}}impl<T>Associative for MultiplicativeOperation<T>where T:Clone+One+Mul<Output=T>{}impl<T>Commutative for MultiplicativeOperation<T>where T:Clone+One+Mul<Output=T>{}impl<T>Invertible for MultiplicativeOperation<T>where T:Clone+One+Mul<Output=T>+Div<Output=T>{#[inline]fn inverse(x:&Self::T)->Self::T{Self::unit().div(x.clone())}#[inline]fn rinv_operate(x:&Self::T,y:&Self::T)->Self::T{(x.clone()).div(y.clone())}}}
// codesnip-guard: algebra
#[cfg_attr(any(),rustfmt::skip)]pub use self::magma::*;#[cfg_attr(any(),rustfmt::skip)]mod magma{#![doc=" algebraic traits"]#[doc=" binary operaion: $T \\circ T \\to T$"]pub trait Magma{#[doc=" type of operands: $T$"]type T:Clone;#[doc=" binary operaion: $\\circ$"]fn operate(x:&Self::T,y:&Self::T)->Self::T;fn reverse_operate(x:&Self::T,y:&Self::T)->Self::T{Self::operate(y,x)}fn operate_assign(x:&mut Self::T,y:&Self::T){*x=Self::operate(x,y);}}#[doc=" $\\forall a,\\forall b,\\forall c \\in T, (a \\circ b) \\circ c = a \\circ (b \\circ c)$"]pub trait Associative:Magma{}#[doc=" associative binary operation"]pub trait SemiGroup:Magma+Associative{}impl<S>SemiGroup for S where S:Magma+Associative{}#[doc=" $\\exists e \\in T, \\forall a \\in T, e \\circ a = a \\circ e = a$"]pub trait Unital:Magma{#[doc=" identity element: $e$"]fn unit()->Self::T;fn is_unit(x:&Self::T)->bool where Self::T:PartialEq{x==&Self::unit()}fn set_unit(x:&mut Self::T){*x=Self::unit();}}pub trait ExpBits{type Iter:Iterator<Item=bool>;fn bits(self)->Self::Iter;}pub trait SignedExpBits{type T:ExpBits;fn neg_and_bits(self)->(bool,Self::T);}pub struct Bits<T>{n:T}macro_rules!impl_exp_bits_for_uint{($($t:ty)*)=>{$(impl Iterator for Bits<$t>{type Item=bool;fn next(&mut self)->Option<bool>{if self.n==0{None}else{let bit=(self.n&1)==1;self.n>>=1;Some(bit)}}}impl ExpBits for$t{type Iter=Bits<$t>;fn bits(self)->Self::Iter{Bits{n:self}}}impl SignedExpBits for$t{type T=$t;fn neg_and_bits(self)->(bool,Self::T){(false,self)}})*};}impl_exp_bits_for_uint!(u8 u16 u32 u64 u128 usize);macro_rules!impl_signed_exp_bits_for_sint{($($s:ty,$u:ty;)*)=>{$(impl SignedExpBits for$s{type T=$u;fn neg_and_bits(self)->(bool,Self::T){(self<0,self.unsigned_abs())}})*};}impl_signed_exp_bits_for_sint!(i8,u8;i16,u16;i32,u32;i64,u64;i128,u128;isize,usize;);#[doc=" associative binary operation and an identity element"]pub trait Monoid:SemiGroup+Unital{#[doc=" binary exponentiation: $x^n = x\\circ\\ddots\\circ x$"]fn pow<E>(mut x:Self::T,exp:E)->Self::T where E:ExpBits{let mut bits=exp.bits();let Some(first)=bits.next()else{return Self::unit();};let mut res=if first{x.clone()}else{Self::unit()};for bit in bits{x=Self::operate(&x,&x);if bit{res=Self::operate(&res,&x);}}res}fn fold<I>(iter:I)->Self::T where I:IntoIterator<Item=Self::T>{let mut iter=iter.into_iter();if let Some(item)=iter.next(){iter.fold(item,|acc,x|Self::operate(&acc,&x))}else{Self::unit()}}}impl<M>Monoid for M where M:SemiGroup+Unital{}#[doc=" $\\exists e \\in T, \\forall a \\in T, \\exists b,c \\in T, b \\circ a = a \\circ c = e$"]pub trait Invertible:Magma+Unital{#[doc=" $a$ where $a \\circ x = e$"]fn inverse(x:&Self::T)->Self::T;fn rinv_operate(x:&Self::T,y:&Self::T)->Self::T{Self::operate(x,&Self::inverse(y))}fn rinv_operate_assign(x:&mut Self::T,y:&Self::T){*x=Self::rinv_operate(x,y);}}#[doc=" associative binary operation and an identity element and inverse elements"]pub trait Group:Monoid+Invertible{fn signed_pow<E>(x:Self::T,exp:E)->Self::T where E:SignedExpBits{let(neg,exp)=E::neg_and_bits(exp);let res=Self::pow(x,exp);if neg{Self::inverse(&res)}else{res}}}impl<G>Group for G where G:Monoid+Invertible{}#[doc=" $\\forall a,\\forall b \\in T, a \\circ b = b \\circ a$"]pub trait Commutative:Magma{}#[doc=" commutative monoid"]pub trait AbelianMonoid:Monoid+Commutative{}impl<M>AbelianMonoid for M where M:Monoid+Commutative{}#[doc=" commutative group"]pub trait AbelianGroup:Group+Commutative{}impl<G>AbelianGroup for G where G:Group+Commutative{}#[doc=" $\\forall a \\in T, a \\circ a = a$"]pub trait Idempotent:Magma{}#[doc=" idempotent monoid"]pub trait IdempotentMonoid:Monoid+Idempotent{}impl<M>IdempotentMonoid for M where M:Monoid+Idempotent{}#[macro_export]macro_rules!monoid_fold{($m:ty)=>{<$m as Unital>::unit()};($m:ty,)=>{<$m as Unital>::unit()};($m:ty,$f:expr)=>{$f};($m:ty,$f:expr,$($ff:expr),*)=>{<$m as Magma>::operate(&($f),&monoid_fold!($m,$($ff),*))};}#[macro_export]macro_rules!define_monoid{($Name:ident,$t:ty,|$x:ident,$y:ident|$op:expr,$unit:expr)=>{struct$Name;impl Magma for$Name{type T=$t;fn operate($x:&Self::T,$y:&Self::T)->Self::T{$op}}impl Unital for$Name{fn unit()->Self::T{$unit}}impl Associative for$Name{}};}}
// codesnip-guard: ring
#[cfg_attr(any(),rustfmt::skip)]pub use self::ring::*;#[cfg_attr(any(),rustfmt::skip)]mod ring{use super::*;use std::{marker::PhantomData,num::Wrapping,ops::{Add,Mul}};pub trait SemiRing{type T:Clone;type Additive:AbelianMonoid<T=Self::T>;type Multiplicative:Monoid<T=Self::T>;#[doc=" additive identity: $0$"]fn zero()->Self::T{<Self::Additive as Unital>::unit()}fn is_zero(x:&Self::T)->bool where Self::T:PartialEq{*x==Self::zero()}#[doc=" multiplicative identity: $1$"]fn one()->Self::T{<Self::Multiplicative as Unital>::unit()}fn is_one(x:&Self::T)->bool where Self::T:PartialEq{*x==Self::one()}#[doc=" additive operaion: $+$"]fn add(x:&Self::T,y:&Self::T)->Self::T{<Self::Additive as Magma>::operate(x,y)}#[doc=" multiplicative operaion: $+$"]fn mul(x:&Self::T,y:&Self::T)->Self::T{<Self::Multiplicative as Magma>::operate(x,y)}fn try_matrix_product(_a:&[Vec<Self::T>],_b:&[Vec<Self::T>])->Option<Vec<Vec<Self::T>>>{None}fn dot_product(x:&[Self::T],y:&[Self::T])->Self::T{assert_eq!(x.len(),y.len());x.iter().zip(y).fold(Self::zero(),|mut sum,(x,y)|{Self::add_assign(&mut sum,&Self::mul(x,y));sum})}fn add_scaled_assign(x:&mut[Self::T],y:&[Self::T],a:&Self::T){assert_eq!(x.len(),y.len());for(x,y)in x.iter_mut().zip(y){Self::add_assign(x,&Self::mul(a,y));}}fn add_assign(x:&mut Self::T,y:&Self::T){<Self::Additive as Magma>::operate_assign(x,y);}fn mul_assign(x:&mut Self::T,y:&Self::T){<Self::Multiplicative as Magma>::operate_assign(x,y);}}pub trait Ring:SemiRing<Additive:Invertible>{#[doc=" additive inverse: $-$"]fn neg(x:&Self::T)->Self::T{<Self::Additive as Invertible>::inverse(x)}#[doc=" additive right inversed operaion: $-$"]fn sub(x:&Self::T,y:&Self::T)->Self::T{<Self::Additive as Invertible>::rinv_operate(x,y)}fn sub_assign(x:&mut Self::T,y:&Self::T){<Self::Additive as Invertible>::rinv_operate_assign(x,y);}}impl<R>Ring for R where R:SemiRing<Additive:Invertible>{}pub trait Field:Ring<Multiplicative:Invertible>{#[doc=" multiplicative inverse: $-$"]fn inv(x:&Self::T)->Self::T{<Self::Multiplicative as Invertible>::inverse(x)}#[doc=" multiplicative right inversed operaion: $-$"]fn div(x:&Self::T,y:&Self::T)->Self::T{<Self::Multiplicative as Invertible>::rinv_operate(x,y)}fn div_assign(x:&mut Self::T,y:&Self::T){<Self::Multiplicative as Invertible>::rinv_operate_assign(x,y);}}impl<F>Field for F where F:Ring<Multiplicative:Invertible>{}pub trait DotProduct:Sized+Clone+Zero+Add<Output=Self>+Mul<Output=Self>{fn try_matrix_product(_a:&[Vec<Self>],_b:&[Vec<Self>])->Option<Vec<Vec<Self>>>{None}fn add_scaled_assign(x:&mut[Self],y:&[Self],a:&Self){assert_eq!(x.len(),y.len());for(x,y)in x.iter_mut().zip(y){*x=x.clone()+a.clone()*y.clone();}}fn dot_product(x:&[Self],y:&[Self])->Self{assert_eq!(x.len(),y.len());x.iter().zip(y).fold(Self::zero(),|sum,(x,y)|sum+x.clone()*y.clone())}}macro_rules!impl_dot_product{($($t:ty)*)=>{$(impl DotProduct for$t{})*};}impl_dot_product!(u8 u16 u32 u64 usize u128 i8 i16 i32 i64 isize i128 f32 f64);impl<T>DotProduct for Wrapping<T>where Wrapping<T>:Clone+Zero+Add<Output=Wrapping<T>>+Mul<Output=Wrapping<T>>{}#[doc=" $+,\\times$"]pub struct AddMulOperation<T>where T:DotProduct+One{_marker:PhantomData<fn()->T>}impl<T>SemiRing for AddMulOperation<T>where T:DotProduct+One{type T=T;type Additive=AdditiveOperation<T>;type Multiplicative=MultiplicativeOperation<T>;#[inline]fn try_matrix_product(a:&[Vec<T>],b:&[Vec<T>])->Option<Vec<Vec<T>>>{T::try_matrix_product(a,b)}#[inline]fn dot_product(x:&[Self::T],y:&[Self::T])->Self::T{T::dot_product(x,y)}#[inline]fn add_scaled_assign(x:&mut[Self::T],y:&[Self::T],a:&Self::T){T::add_scaled_assign(x,y,a);}}}
// codesnip-guard: zero_one
#[cfg_attr(any(),rustfmt::skip)]pub use self::zero_one::{One,Zero};#[cfg_attr(any(),rustfmt::skip)]mod zero_one{pub trait Zero:Sized{fn zero()->Self;#[inline]fn is_zero(&self)->bool where Self:PartialEq{self==&Self::zero()}#[inline]fn set_zero(&mut self){*self=Self::zero();}}pub trait One:Sized{fn one()->Self;#[inline]fn is_one(&self)->bool where Self:PartialEq{self==&Self::one()}#[inline]fn set_one(&mut self){*self=Self::one();}}macro_rules!impl_zero_one{($({$(<$T:ident:$Bound:ident>)?$Trait:ident$method:ident$($t:ty)*,$e:expr})*)=>{$(impl_zero_one!(@impl[$(<$T:$Bound>)?]$Trait$method[$($t)*],$e);)*};(@impl[<$T:ident:$Bound:ident>]$Trait:ident$method:ident[$($t:ty)*],$e:expr)=>{$(impl<$T:$Bound>$Trait for$t{fn$method()->Self{$e}})*};(@impl[]$Trait:ident$method:ident[$($t:ty)*],$e:expr)=>{$(impl$Trait for$t{fn$method()->Self{$e}})*};}impl_zero_one!({Zero zero u8 u16 u32 u64 usize i8 i16 i32 i64 isize u128 i128,0}{Zero zero f32 f64,0.}{One one u8 u16 u32 u64 usize i8 i16 i32 i64 isize u128 i128,1}{One one f32 f64,1.}{<T:Zero>Zero zero std::num::Wrapping<T>,Self(T::zero())}{<T:One>One one std::num::Wrapping<T>,Self(T::one())});}