// codesnip-guard: automata_learning #[cfg_attr(any(),rustfmt::skip)]pub use self::automata_learning::*;#[cfg_attr(any(),rustfmt::skip)]mod automata_learning{use super::{BitSet,Field,Invertible,Matrix,RandomSpec,SerdeByteStr,Xorshift};use std::{cell::RefCell,collections::{HashMap,HashSet},fmt::{self,Debug},iter::{from_fn,once_with},marker::PhantomData,time::Instant};pub trait BlackBoxAutomaton{type Output;fn sigma(&self)->usize;fn behavior(&self,input:I)->Self::Output where I:IntoIterator;}#[derive(Debug,Clone)]pub struct BlackBoxAutomatonImplwhere F:Fn(Vec)->T{sigma:usize,behavior_fn:F,memo:RefCell,T>>}implBlackBoxAutomatonImplwhere F:Fn(Vec)->T{pub fn new(sigma:usize,behavior_fn:F)->Self{Self{sigma,behavior_fn,memo:RefCell::new(HashMap::new())}}}implBlackBoxAutomaton for BlackBoxAutomatonImplwhere F:Fn(Vec)->T,T:Clone{type Output=T;fn sigma(&self)->usize{self.sigma}fn behavior(&self,input:I)->Self::Output where I:IntoIterator{let input:Vec=input.into_iter().collect();self.memo.borrow_mut().entry(input.clone()).or_insert_with(| |(self.behavior_fn)(input)).clone()}}implBlackBoxAutomaton for&A where A:BlackBoxAutomaton{type Output=A::Output;fn sigma(&self)->usize{(*self).sigma()}fn behavior(&self,input:I)->Self::Output where I:IntoIterator{(*self).behavior(input)}}#[derive(Debug,Clone)]struct DfaState{delta:Vec,accept:bool}#[derive(Debug,Clone)]pub struct DeterministicFiniteAutomaton{states:Vec,initial_state:usize}impl DeterministicFiniteAutomaton{pub fn size(&self)->usize{self.states.len()}pub fn delta(&self,state:usize,input:usize)->usize{assert!(statebool{assert!(stateusize{self.states[0].delta.len()}fn behavior(&self,input:I)->Self::Output where I:IntoIterator{let mut state=self.initial_state;for x in input{state=self.states[state].delta[x];}self.states[state].accept}}impl SerdeByteStr for DfaState{fn serialize(&self,buf:&mut Vec){self.delta.serialize(buf);self.accept.serialize(buf);}fn deserialize(iter:&mut I)->Self where I:Iterator{let delta=Vec::deserialize(iter);let accept=bool::deserialize(iter);Self{delta,accept}}}impl SerdeByteStr for DeterministicFiniteAutomaton{fn serialize(&self,buf:&mut Vec){self.states.serialize(buf);self.initial_state.serialize(buf);}fn deserialize(iter:&mut I)->Self where I:Iterator{let states=Vec::deserialize(iter);let initial_state=usize::deserialize(iter);Self{states,initial_state}}}pub struct WeightedFiniteAutomatonwhere F:Field{pub initial_weights:Matrix,pub transitions:Vec>,pub final_weights:Matrix}implDebug for WeightedFiniteAutomatonwhere F:Field{fn fmt(&self,f:&mut fmt::Formatter<'_>)->fmt::Result{f.debug_struct("WeightedFiniteAutomaton").field("initial_weights",&self.initial_weights).field("transitions",&self.transitions).field("final_weights",&self.final_weights).finish()}}implClone for WeightedFiniteAutomatonwhere F:Field{fn clone(&self)->Self{Self{initial_weights:self.initial_weights.clone(),transitions:self.transitions.clone(),final_weights:self.final_weights.clone()}}}implBlackBoxAutomaton for WeightedFiniteAutomatonwhere F:Field{type Output=F::T;fn sigma(&self)->usize{self.transitions.len()}fn behavior(&self,input:I)->Self::Output where I:IntoIterator{let mut weights=self.initial_weights.clone();for x in input{weights=&weights*&self.transitions[x];}let result=&weights*&self.final_weights;if result.shape!=(0,0){result[0][0].clone()}else{F::zero()}}}implSerdeByteStr for WeightedFiniteAutomatonwhere F:Field{fn serialize(&self,buf:&mut Vec){self.initial_weights.serialize(buf);self.transitions.serialize(buf);self.final_weights.serialize(buf);}fn deserialize(iter:&mut I)->Self where I:Iterator{let initial_weights=Matrix::deserialize(iter);let transitions=Vec::deserialize(iter);let final_weights=Matrix::deserialize(iter);Self{initial_weights,transitions,final_weights}}}pub fn dense_sampling(sigma:usize,max_len:usize)->impl Iterator>{assert_ne!(sigma,0,"Sigma must be greater than 0");let mut current=vec![];once_with(Vec::new).chain(from_fn(move| |{let mut carry=true;for i in(0..current.len()).rev(){current[i]+=1;if current[i]==sigma{current[i]=0;}else{carry=false;break;}}if carry{current.push(0);}if current.len()<=max_len{Some(current.to_vec())}else{None}}))}pub fn random_sampling(sigma:usize,len_spec:impl RandomSpec,seconds:f64)->impl Iterator>{assert_ne!(sigma,0,"Sigma must be greater than 0");let now=Instant::now();let mut rng=Xorshift::new();from_fn(move| |{if now.elapsed().as_secs_f64()>seconds{None}else{let n=rng.random(&len_spec);Some(rng.random_iter(0..sigma).take(n).collect())}})}#[derive(Debug,Clone)]pub struct DfaLearningwhere A:BlackBoxAutomaton{automaton:A,prefixes:Vec>,suffixes:Vec>,table:Vec,row_map:HashMap}implDfaLearningwhere A:BlackBoxAutomaton{pub fn new(automaton:A)->Self{let mut this=Self{automaton,prefixes:vec![],suffixes:vec![],table:vec![],row_map:HashMap::new()};this.add_suffix(vec![]);this.add_prefix(vec![]);this}fn add_prefix(&mut self,prefix:Vec)->usize{let row:BitSet=self.suffixes.iter().map(|s|{self.automaton.behavior(prefix.iter().cloned().chain(s.iter().cloned()))}).collect();*self.row_map.entry(row.clone()).or_insert_with(| |{let idx=self.table.len();self.table.push(row);self.prefixes.push(prefix);idx})}fn add_suffix(&mut self,suffix:Vec){if self.suffixes.contains(&suffix){return;}for(prefix,table)in self.prefixes.iter_mut().zip(&mut self.table){table.push(self.automaton.behavior(prefix.iter().cloned().chain(suffix.iter().cloned())));}self.suffixes.push(suffix);self.row_map.clear();for(i_prefix,row)in self.table.iter().enumerate(){self.row_map.insert(row.clone(),i_prefix);}}pub fn construct_dfa(&mut self)->DeterministicFiniteAutomaton{let sigma=self.automaton.sigma();let mut dfa=DeterministicFiniteAutomaton{states:vec![],initial_state:0};let mut i_prefix=0;while i_prefix=self.prefixes[i_prefix].iter().cloned().chain([x]).collect();let index=self.add_prefix(prefix);delta.push(index);}dfa.states.push(DfaState{delta,accept:self.table[i_prefix].get(0)});i_prefix+=1;}dfa}pub fn train_sample(&mut self,dfa:&DeterministicFiniteAutomaton,sample:&[usize])->bool{let expected=self.automaton.behavior(sample.iter().cloned());if expected==dfa.behavior(sample.iter().cloned()){return false;}let n=sample.len();let mut states:Vec<(usize,usize)>=Vec::with_capacity(n+1);let mut s=0usize;states.push((s,0));for(k,&x)in sample.iter().enumerate(){s=dfa.states[s].delta[x];states.push((s,k+1));}let split=states.partition_point(|&(state,k)|{self.automaton.behavior(self.prefixes[state].iter().cloned().chain(sample[k..].iter().cloned()))==expected});let new_prefix=sample[..split].to_vec();let new_suffix=sample[split..].to_vec();self.add_suffix(new_suffix);self.add_prefix(new_prefix);true}pub fn train(&mut self,samples:impl IntoIterator>)->DeterministicFiniteAutomaton{let mut dfa=self.construct_dfa();for sample in samples{if self.train_sample(&dfa,&sample){dfa=self.construct_dfa();}}dfa}}pub struct WfaLearningwhere F:Field,A:BlackBoxAutomaton{automaton:A,prefixes:Vec>,suffixes:Vec>,inv_h:Matrix,nh:Vec>,wfa:WeightedFiniteAutomaton,_marker:PhantomDataF>}implDebug for WfaLearningwhere F:Field,A:BlackBoxAutomaton+Debug{fn fmt(&self,f:&mut fmt::Formatter<'_>)->fmt::Result{f.debug_struct("WfaLearning").field("automaton",&self.automaton).field("prefixes",&self.prefixes).field("suffixes",&self.suffixes).field("inv_h",&self.inv_h).field("nh",&self.nh).field("wfa",&self.wfa).finish()}}implClone for WfaLearningwhere F:Field,A:BlackBoxAutomaton+Clone{fn clone(&self)->Self{Self{automaton:self.automaton.clone(),prefixes:self.prefixes.clone(),suffixes:self.suffixes.clone(),inv_h:self.inv_h.clone(),nh:self.nh.clone(),wfa:self.wfa.clone(),_marker:self._marker}}}implWfaLearningwhere F:Field,A:BlackBoxAutomaton{pub fn new(automaton:A)->Self{let sigma=automaton.sigma();Self{automaton,prefixes:vec![],suffixes:vec![],inv_h:Matrix::zeros((0,0)),nh:vec![Matrix::zeros((0,0));sigma],wfa:WeightedFiniteAutomaton{initial_weights:Matrix::zeros((1,0)),transitions:vec![Matrix::zeros((0,0));sigma],final_weights:Matrix::zeros((0,1))},_marker:PhantomData}}pub fn wfa(&self)->&WeightedFiniteAutomaton{&self.wfa}fn split_sample(&mut self,sample:&[usize])->Option<(Vec,Vec)>{if self.prefixes.is_empty()&&!F::is_zero(&self.automaton.behavior(sample.iter().cloned())){return Some((vec![],sample.to_vec()));}let expected=self.automaton.behavior(sample.iter().cloned());if expected==self.wfa.behavior(sample.iter().cloned()){return None;}let n=sample.len();let dim=self.wfa.final_weights.shape.0;let mut states:Vec<(Matrix,usize)>=Vec::with_capacity(n+1);let mut v=self.wfa.final_weights.clone();states.push((v.clone(),n));for k in(0..n).rev(){v=&self.wfa.transitions[sample[k]]*&v;states.push((v.clone(),k));}states.reverse();let split=states.partition_point(|(state,k)|{(0..dim).any(|j|{self.automaton.behavior(self.prefixes[j].iter().cloned().chain(sample[*k..].iter().cloned()))!=state[j][0]})});Some((sample[..split].to_vec(),sample[split..].to_vec()))}pub fn train_sample(&mut self,sample:&[usize])->bool{let Some((prefix,suffix))=self.split_sample(sample)else{return false;};self.prefixes.push(prefix);self.suffixes.push(suffix);let n=self.inv_h.shape.0;let prefix=&self.prefixes[n];let suffix=&self.suffixes[n];let u=Matrix::::new_with((n,1),|i,_|{self.automaton.behavior(self.prefixes[i].iter().cloned().chain(suffix.iter().cloned()))});let v=Matrix::::new_with((1,n),|_,j|{self.automaton.behavior(prefix.iter().cloned().chain(self.suffixes[j].iter().cloned()))});let w=Matrix::::new_with((1,1),|_,_|{self.automaton.behavior(prefix.iter().cloned().chain(suffix.iter().cloned()))});let t=&self.inv_h*&u;let s=&v*&self.inv_h;let d=F::inv(&(&w-&(&v*&t))[0][0]);let dh=&t*&s;for i in 0..n{for j in 0..n{F::add_assign(&mut self.inv_h[i][j],&F::mul(&dh[i][j],&d));}}self.inv_h.add_col_with(|i,_|F::neg(&F::mul(&t[i][0],&d)));self.inv_h.add_row_with(|_,j|{if j!=n{F::neg(&F::mul(&s[0][j],&d))}else{d.clone()}});for(x,transition)in self.wfa.transitions.iter_mut().enumerate(){let b=&(&self.nh[x]*&t)*&s;for i in 0..n{for j in 0..n{F::add_assign(&mut transition[i][j],&F::mul(&b[i][j],&d));}}}for(x,nh)in self.nh.iter_mut().enumerate(){nh.add_col_with(|i,j|{self.automaton.behavior(self.prefixes[i].iter().cloned().chain([x]).chain(self.suffixes[j].iter().cloned()))});nh.add_row_with(|i,j|{self.automaton.behavior(self.prefixes[i].iter().cloned().chain([x]).chain(self.suffixes[j].iter().cloned()))});}self.wfa.initial_weights.add_col_with(|_,_|if n==0{F::one()}else{F::zero()});self.wfa.final_weights.add_row_with(|_,_|self.automaton.behavior(prefix.iter().cloned()));for(x,transition)in self.wfa.transitions.iter_mut().enumerate(){transition.add_col_with(|_,_|F::zero());transition.add_row_with(|_,_|F::zero());for i in 0..=n{for j in 0..=n{if i==n||j==n{for k in 0..=n{if i!=n&&j!=n&&k!=n{continue;}F::add_assign(&mut transition[i][k],&F::mul(&self.nh[x][i][j],&self.inv_h[j][k]));}}else{let k=n;F::add_assign(&mut transition[i][k],&F::mul(&self.nh[x][i][j],&self.inv_h[j][k]));}}}}true}pub fn train(&mut self,samples:impl IntoIterator>){for sample in samples{self.train_sample(&sample);}}pub fn batch_train(&mut self,samples:impl IntoIterator>){let mut prefix_set:HashSet<_>=self.prefixes.iter().cloned().collect();let mut suffix_set:HashSet<_>=self.suffixes.iter().cloned().collect();for sample in samples{if prefix_set.insert(sample.to_vec()){self.prefixes.push(sample.to_vec());}if suffix_set.insert(sample.to_vec()){self.suffixes.push(sample);}}let mut h=Matrix::::new_with((self.prefixes.len(),self.suffixes.len()),|i,j|{self.automaton.behavior(self.prefixes[i].iter().cloned().chain(self.suffixes[j].iter().cloned()))});if!self.prefixes.is_empty()&&!self.suffixes.is_empty()&&F::is_zero(&h[0][0]){for j in 1..self.suffixes.len(){if!F::is_zero(&h[0][j]){self.suffixes.swap(0,j);for row in&mut h.data{row.swap(0,j);}break;}}}let mut row_id:Vec=(0..h.shape.0).collect();let mut pivots=vec![];h.row_reduction_with(false,|r,p,c|{row_id.swap(r,p);pivots.push((row_id[r],c));});let mut new_prefixes=vec![];let mut new_suffixes=vec![];for(i,j)in pivots{new_prefixes.push(self.prefixes[i].clone());new_suffixes.push(self.suffixes[j].clone());}self.prefixes=new_prefixes;self.suffixes=new_suffixes;assert_eq!(self.prefixes.len(),self.suffixes.len());let n=self.prefixes.len();let h=Matrix::::new_with((n,n),|i,j|{self.automaton.behavior(self.prefixes[i].iter().cloned().chain(self.suffixes[j].iter().cloned()))});self.inv_h=h.inverse().expect("Hankel matrix must be invertible");self.wfa=WeightedFiniteAutomaton::{initial_weights:Matrix::new_with((1,n),|_,j|{if self.prefixes[j].is_empty(){F::one()}else{F::zero()}}),transitions:(0..self.automaton.sigma()).map(|x|{&Matrix::new_with((n,n),|i,j|{self.automaton.behavior(self.prefixes[i].iter().cloned().chain([x]).chain(self.suffixes[j].iter().cloned()))})*&self.inv_h}).collect(),final_weights:Matrix::new_with((n,1),|i,_|{self.automaton.behavior(self.prefixes[i].iter().cloned())})};}}} // 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 AdditiveOperationwhere T:Clone+Zero+Add{_marker:PhantomDataT>}implMagma for AdditiveOperationwhere T:Clone+Zero+Add{type T=T;#[inline]fn operate(x:&Self::T,y:&Self::T)->Self::T{x.clone()+y.clone()}}implUnital for AdditiveOperationwhere T:Clone+Zero+Add{#[inline]fn unit()->Self::T{Zero::zero()}}implAssociative for AdditiveOperationwhere T:Clone+Zero+Add{}implCommutative for AdditiveOperationwhere T:Clone+Zero+Add{}implInvertible for AdditiveOperationwhere T:Clone+Zero+Add+Sub+Neg{#[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: BitSet #[cfg_attr(any(),rustfmt::skip)]pub use self::bitset::BitSet;#[cfg_attr(any(),rustfmt::skip)]mod bitset{#[cfg(target_arch="x86_64")]use super::avx512_enabled;use std::{cmp::Ordering,ops::{BitAnd,BitAndAssign,BitOr,BitOrAssign,BitXor,BitXorAssign,Not,Shl,ShlAssign,Shr,ShrAssign}};const BIT_AND:u8=0;const BIT_OR:u8=1;const BIT_XOR:u8=2;#[cfg(target_arch="x86_64")]const SIMD_MIN_BLOCKS:usize=8;#[repr(C,align(64))]#[derive(Clone,Copy,Debug,Default,PartialEq,Eq,PartialOrd,Ord,Hash)]struct Block([u64;8]);#[derive(Clone,Debug,Default,PartialEq,Eq,PartialOrd,Ord,Hash)]pub struct BitSet{size:usize,bits:Vec}impl BitSet{pub fn new(size:usize)->Self{Self{size,bits:vec![Block::default();size.div_ceil(512)]}}pub fn len(&self)->usize{self.size}pub fn is_empty(&self)->bool{self.size==0}#[doc=" Parses ASCII `0` and `1`, with the first character at bit index zero."]#[doc=" Returns `None` if any other character occurs."]pub fn from_binary(s:&str)->Option{let bytes=s.as_bytes();let mut bits=Self::new(bytes.len());let end=bytes.len()/64*64;#[cfg(target_arch="x86_64")]let parsed=if avx512_enabled()&&is_x86_feature_detected!("avx512bw"){unsafe{simd::parse_binary_avx512(&bytes[..end],bits.words_mut())}}else if is_x86_feature_detected!("avx2"){unsafe{simd::parse_binary_avx2(&bytes[..end],bits.words_mut())}}else{Self::parse_binary_scalar(&bytes[..end],bits.words_mut())};#[cfg(not(target_arch="x86_64"))]let parsed=Self::parse_binary_scalar(&bytes[..end],bits.words_mut());if!parsed{return None;}if end!=bytes.len(){let mut word=0;for(i,&b)in bytes[end..].iter().enumerate(){if b!=b'0'&&b!=b'1'{return None;}word|=u64::from(b&1)<bool{for(chunk,word)in bytes.as_chunks::<64>().0 .iter().zip(words){for(i,byte)in chunk.as_chunks::<8>().0 .iter().enumerate(){let x=u64::from_le_bytes(*byte);if x&0xfefe_fefe_fefe_fefe!=0x3030_3030_3030_3030{return false;}*word|=((x&0x0101_0101_0101_0101).wrapping_mul(0x0102_0408_1020_4080)>>56)<<(i*8);}}true}#[doc=" Returns ASCII `0` and `1` in increasing bit-index order."]pub fn to_binary(&self)->String{const TABLE:[[u8;8];256]={let mut table=[[b'0';8];256];let mut i=0;while i<256{let mut j=0;while j<8{table[i][j]|=((i>>j)&1)as u8;j+=1;}i+=1;}table};let mut bytes=vec![b'0';self.size.div_ceil(8)*8];#[cfg(target_arch="x86_64")]let end=if self.size>=64&&avx512_enabled()&&is_x86_feature_detected!("avx512bw"){let end=self.size/64*64;unsafe{simd::write_binary_avx512(&mut bytes[..end],self.words());}end}else if self.size>=64&&is_x86_feature_detected!("avx2"){let end=self.size/64*64;unsafe{simd::write_binary_avx2(&mut bytes[..end],self.words());}end}else{0};#[cfg(not(target_arch="x86_64"))]let end=0;for(chunk,&word)in bytes[end..].chunks_mut(64).zip(&self.words()[end/64..]){for(i,byte)in chunk.as_chunks_mut::<8>().0 .iter_mut().enumerate(){byte.copy_from_slice(&TABLE[(word>>(i*8)&255)as usize]);}}bytes.truncate(self.size);unsafe{String::from_utf8_unchecked(bytes)}}pub fn ones(size:usize)->Self{let mut self_=Self{size,bits:vec![Block([u64::MAX;8]);size.div_ceil(512)]};self_.trim();self_}pub fn get(&self,i:usize)->bool{self.bits[i>>9].0[i>>6&7]&(1<<(i&63))!=0}pub fn set(&mut self,i:usize,b:bool){let word=&mut self.bits[i>>9].0[i>>6&7];if b{*word|=1<<(i&63);}else{*word&=!(1<<(i&63));}}#[doc=" Clears all bits."]pub fn reset(&mut self){self.bits.fill(Block::default());}#[doc=" Sets all bits to `value`."]pub fn fill(&mut self,value:bool){self.bits.fill(Block([if value{u64::MAX}else{0};8]));self.trim();}#[doc=" Tests whether any bit is set."]#[inline]pub fn any(&self)->bool{!self.none()}#[doc=" Tests whether all bits are unset."]#[inline]pub fn none(&self)->bool{#[cfg(target_arch="x86_64")]if self.bits.len()>=SIMD_MIN_BLOCKS{if self.bits[0].0[0]!=0{return false;}if avx512_enabled()&&is_x86_feature_detected!("avx512f"){return unsafe{simd::none_avx512(&self.bits)};}if is_x86_feature_detected!("avx2"){return unsafe{simd::none_avx2(&self.bits)};}}self.words().iter().all(|&word|word==0)}#[doc=" Tests whether all bits are set."]#[inline(always)]pub fn all(&self)->bool{let words=self.words();let full_words=self.size>>6;#[cfg(target_arch="x86_64")]if self.size>>9>=SIMD_MIN_BLOCKS{let full_blocks=self.size>>9;if words[0]!=u64::MAX{return false;}let full_blocks_are_set=if avx512_enabled()&&is_x86_feature_detected!("avx512f"){Some(unsafe{simd::all_avx512(&self.bits[..full_blocks])})}else if is_x86_feature_detected!("avx2"){Some(unsafe{simd::all_avx2(&self.bits[..full_blocks])})}else{None};if let Some(full_blocks_are_set)=full_blocks_are_set{return full_blocks_are_set&&words[full_blocks*8..full_words].iter().all(|&word|word==u64::MAX)&&(self.size&63==0||words[full_words]==u64::MAX>>(64-(self.size&63)));}}if self.size&63==0{return words.iter().all(|&word|word==u64::MAX);}words[..full_words].iter().all(|&word|word==u64::MAX)&&words[full_words]==u64::MAX>>(64-(self.size&63))}#[doc=" Iterates over set-bit indices in ascending order."]pub fn iter_ones(&self)->impl Iterator+'_{self.words().iter().copied().enumerate().flat_map(|(word_index,mut word)|{std::iter::from_fn(move| |{if word==0{None}else{let bit=word.trailing_zeros()as usize;word&=word-1;Some((word_index<<6)|bit)}})})}#[doc=" Counts set bits."]#[inline]pub fn count_ones(&self)->u64{let words=self.words();#[cfg(target_arch="x86_64")]if words.len()>=8{if avx512_enabled()&&is_x86_feature_detected!("avx512f")&&is_x86_feature_detected!("avx512vpopcntdq"){return unsafe{simd::count_ones_avx512(&self.bits)};}if is_x86_feature_detected!("avx2"){return unsafe{simd::count_ones_avx2(&self.bits)};}}words.iter().map(|word|word.count_ones()as u64).sum()}#[doc=" Counts unset bits."]#[inline]pub fn count_zeros(&self)->u64{self.size as u64-self.count_ones()}pub fn push(&mut self,b:bool){if self.size&511==0{self.bits.push(Block::default());}if b{self.bits[self.size>>9].0[self.size>>6&7]|=1<<(self.size&63);}self.size+=1;}pub fn resize(&mut self,new_size:usize){match self.size.cmp(&new_size){Ordering::Less=>self.bits.resize(new_size.div_ceil(512),Block::default()),Ordering::Equal=>{}Ordering::Greater=>self.bits.truncate(new_size.div_ceil(512)),}self.size=new_size;self.trim();}#[doc=" Assigns `self | (self << rhs)` to `self`."]#[inline]pub fn shl_bitor_assign(&mut self,rhs:usize){self.shift_left::(rhs);}#[doc=" Assigns `self | (self >> rhs)` to `self`."]#[inline]pub fn shr_bitor_assign(&mut self,rhs:usize){self.shift_right::(rhs);}#[doc=" Returns words in increasing bit-index order; unused high bits in the last word are zero."]pub fn words(&self)->&[u64]{unsafe{std::slice::from_raw_parts(self.bits.as_ptr().cast(),self.size.div_ceil(64))}}#[doc=" Returns mutable words in increasing bit-index order."]#[doc=" Callers must keep unused high bits in the last word zero."]pub fn words_mut(&mut self)->&mut[u64]{let len=self.size.div_ceil(64);unsafe{std::slice::from_raw_parts_mut(self.bits.as_mut_ptr().cast(),len)}}fn trim(&mut self){let used_words=self.size.div_ceil(64)&7;if let Some(last)=self.bits.last_mut(){if self.size&63!=0{last.0[(self.size>>6)&7]&=u64::MAX>>(64-(self.size&63));}if used_words!=0{last.0[used_words..].fill(0);}}}#[inline]fn bitop_assign(&mut self,rhs:&Self){assert_eq!(self.size,rhs.size,"bitwise operations require equal lengths");#[cfg(target_arch="x86_64")]if self.bits.len()>=SIMD_MIN_BLOCKS{if avx512_enabled()&&is_x86_feature_detected!("avx512f"){unsafe{simd::bitop_avx512::(&mut self.bits,&rhs.bits)};return;}if is_x86_feature_detected!("avx2"){unsafe{simd::bitop_avx2::(&mut self.bits,&rhs.bits)};return;}}for(lhs,&rhs)in self.words_mut().iter_mut().zip(rhs.words()){*lhs=match OP{BIT_AND=>*lhs&rhs,BIT_OR=>*lhs|rhs,BIT_XOR=>*lhs^rhs,_=>unreachable!(),};}}#[inline]fn shift_left(&mut self,rhs:usize){if rhs==0{return;}if rhs>=self.size{if!OR_ASSIGN{self.reset();}return;}#[cfg(target_arch="x86_64")]if self.bits.len()>=SIMD_MIN_BLOCKS{if avx512_enabled()&&is_x86_feature_detected!("avx512f")&&is_x86_feature_detected!("avx512vbmi2"){unsafe{simd::shift_left_avx512::(&mut self.bits,rhs)};if self.size&511!=0{self.trim();}return;}if is_x86_feature_detected!("avx2"){unsafe{simd::shift_left_avx2::(self.words_mut(),rhs)};if self.size&63!=0{self.trim();}return;}}let bits=self.words_mut();let word_shift=rhs>>6;let bit_shift=rhs&63;if bit_shift==0{for i in(0..bits.len()-word_shift).rev(){if OR_ASSIGN{bits[i+word_shift]|=bits[i];}else{bits[i+word_shift]=bits[i];}}}else{for i in(1..bits.len()-word_shift).rev(){let value=(bits[i]<>(64-bit_shift));if OR_ASSIGN{bits[i+word_shift]|=value;}else{bits[i+word_shift]=value;}}if OR_ASSIGN{bits[word_shift]|=bits[0]<(&mut self,rhs:usize){if rhs==0{return;}if rhs>=self.size{if!OR_ASSIGN{self.reset();}return;}#[cfg(target_arch="x86_64")]if self.bits.len()>=SIMD_MIN_BLOCKS{if avx512_enabled()&&is_x86_feature_detected!("avx512f")&&is_x86_feature_detected!("avx512vbmi2"){unsafe{simd::shift_right_avx512::(&mut self.bits,rhs)};return;}if is_x86_feature_detected!("avx2"){unsafe{simd::shift_right_avx2::(self.words_mut(),rhs)};return;}}let bits=self.words_mut();let word_shift=rhs>>6;let bit_shift=rhs&63;if bit_shift==0{for i in word_shift..bits.len(){if OR_ASSIGN{bits[i-word_shift]|=bits[i];}else{bits[i-word_shift]=bits[i];}}}else{for i in word_shift..bits.len()-1{let value=(bits[i]>>bit_shift)|(bits[i+1]<<(64-bit_shift));if OR_ASSIGN{bits[i-word_shift]|=value;}else{bits[i-word_shift]=value;}}if OR_ASSIGN{bits[bits.len()-word_shift-1]|=bits[bits.len()-1]>>bit_shift;}else{bits[bits.len()-word_shift-1]=bits[bits.len()-1]>>bit_shift;}}if!OR_ASSIGN{let end=bits.len()-word_shift;bits[end..].fill(0);}}}impl Extendfor BitSet{fn extend>(&mut self,iter:T){for bit in iter{self.push(bit);}}}impl FromIteratorfor BitSet{fn from_iter>(iter:T)->Self{let mut set=BitSet::new(0);set.extend(iter);set}}impl ShlAssignfor BitSet{#[inline]fn shl_assign(&mut self,rhs:usize){self.shift_left::(rhs);}}impl Shlfor BitSet{type Output=Self;fn shl(mut self,rhs:usize)->Self::Output{self<<=rhs;self}}impl ShrAssignfor BitSet{#[inline]fn shr_assign(&mut self,rhs:usize){self.shift_right::(rhs);}}impl Shrfor BitSet{type Output=Self;fn shr(mut self,rhs:usize)->Self::Output{self>>=rhs;self}}impl BitOrAssign<&BitSet>for BitSet{#[inline]fn bitor_assign(&mut self,rhs:&Self){self.bitop_assign::(rhs);}}impl BitOr<&BitSet>for BitSet{type Output=Self;fn bitor(mut self,rhs:&Self)->Self::Output{self|=rhs;self}}impl BitOr<&BitSet>for&BitSet{type Output=BitSet;fn bitor(self,rhs:&BitSet)->Self::Output{let mut res=self.clone();res|=rhs;res}}impl BitAndAssign<&BitSet>for BitSet{#[inline]fn bitand_assign(&mut self,rhs:&Self){self.bitop_assign::(rhs);}}impl BitAnd<&BitSet>for BitSet{type Output=Self;fn bitand(mut self,rhs:&Self)->Self::Output{self&=rhs;self}}impl BitAnd<&BitSet>for&BitSet{type Output=BitSet;fn bitand(self,rhs:&BitSet)->Self::Output{let mut res=self.clone();res&=rhs;res}}impl BitXorAssign<&BitSet>for BitSet{#[inline]fn bitxor_assign(&mut self,rhs:&Self){self.bitop_assign::(rhs);}}impl BitXor<&BitSet>for BitSet{type Output=Self;fn bitxor(mut self,rhs:&Self)->Self::Output{self^=rhs;self}}impl BitXor<&BitSet>for&BitSet{type Output=BitSet;fn bitxor(self,rhs:&BitSet)->Self::Output{let mut res=self.clone();res^=rhs;res}}impl Not for BitSet{type Output=Self;fn not(mut self)->Self::Output{for word in self.words_mut(){*word=!*word;}self.trim();self}}impl Not for&BitSet{type Output=BitSet;fn not(self)->Self::Output{!self.clone()}}#[cfg(target_arch="x86_64")]#[allow(unsafe_op_in_unsafe_fn)]mod simd{#[cfg(target_arch="x86_64")]use super::{BIT_AND,BIT_OR,BIT_XOR,Block};#[cfg(target_arch="x86_64")]use std::arch::x86_64::*;#[target_feature(enable="avx512bw")]#[cfg(target_arch="x86_64")]pub unsafe fn write_binary_avx512(bytes:&mut[u8],words:&[u64]){let one=_mm512_set1_epi8(1);let zero=_mm512_set1_epi8(b'0' as i8);for(chunk,&word)in bytes.as_chunks_mut::<64>().0 .iter_mut().zip(words){let value=_mm512_mask_add_epi8(zero,word,zero,one);_mm512_storeu_si512(chunk.as_mut_ptr().cast(),value);}}#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]pub unsafe fn write_binary_avx2(bytes:&mut[u8],words:&[u64]){let indices=_mm256_setr_epi64x(0,0x0101_0101_0101_0101,0x0202_0202_0202_0202,0x0303_0303_0303_0303);let mask=_mm256_set1_epi64x(0x8040_2010_0804_0201u64 as i64);let one=_mm256_set1_epi8(b'1' as i8);let zero=_mm256_setzero_si256();for(chunk,&word)in bytes.as_chunks_mut::<64>().0 .iter_mut().zip(words){let lo=_mm256_shuffle_epi8(_mm256_set1_epi32(word as i32),indices);let hi=_mm256_shuffle_epi8(_mm256_set1_epi32((word>>32)as i32),indices);let lo=_mm256_add_epi8(one,_mm256_cmpeq_epi8(_mm256_and_si256(lo,mask),zero));let hi=_mm256_add_epi8(one,_mm256_cmpeq_epi8(_mm256_and_si256(hi,mask),zero));_mm256_storeu_si256(chunk.as_mut_ptr().cast(),lo);_mm256_storeu_si256(chunk.as_mut_ptr().add(32).cast(),hi);}}#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]pub unsafe fn parse_binary_avx2(bytes:&[u8],words:&mut[u64])->bool{let one=_mm256_set1_epi8(b'1' as i8);let mask=_mm256_set1_epi8(!1);let zero=_mm256_set1_epi8(b'0' as i8);for(chunk,word)in bytes.as_chunks::<64>().0 .iter().zip(words){let a=_mm256_loadu_si256(chunk.as_ptr().cast());let b=_mm256_loadu_si256(chunk.as_ptr().add(32).cast());let valid_a=_mm256_cmpeq_epi8(_mm256_and_si256(a,mask),zero);let valid_b=_mm256_cmpeq_epi8(_mm256_and_si256(b,mask),zero);if _mm256_movemask_epi8(_mm256_and_si256(valid_a,valid_b))!=-1{return false;}*word=_mm256_movemask_epi8(_mm256_cmpeq_epi8(a,one))as u32 as u64|((_mm256_movemask_epi8(_mm256_cmpeq_epi8(b,one))as u32 as u64)<<32);}true}#[target_feature(enable="avx512bw")]#[cfg(target_arch="x86_64")]pub unsafe fn parse_binary_avx512(bytes:&[u8],words:&mut[u64])->bool{let one=_mm512_set1_epi8(b'1' as i8);let mask=_mm512_set1_epi8(!1);let zero=_mm512_set1_epi8(b'0' as i8);for(chunk,word)in bytes.as_chunks::<64>().0 .iter().zip(words){let x=_mm512_loadu_si512(chunk.as_ptr().cast());if _mm512_cmpeq_epi8_mask(_mm512_and_si512(x,mask),zero)!=u64::MAX{return false;}*word=_mm512_cmpeq_epi8_mask(x,one);}true}#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]pub unsafe fn bitop_avx2(lhs:&mut[Block],rhs:&[Block]){let lhs_ptr=lhs.as_mut_ptr().cast::<__m256i>();let rhs_ptr=rhs.as_ptr().cast::<__m256i>();for i in 0..lhs.len()*2{let lhs_value=_mm256_load_si256(lhs_ptr.add(i));let rhs_value=_mm256_load_si256(rhs_ptr.add(i));let value=match OP{BIT_AND=>_mm256_and_si256(lhs_value,rhs_value),BIT_OR=>_mm256_or_si256(lhs_value,rhs_value),BIT_XOR=>_mm256_xor_si256(lhs_value,rhs_value),_=>unreachable!(),};_mm256_store_si256(lhs_ptr.add(i),value);}}#[target_feature(enable="avx512f")]#[cfg(target_arch="x86_64")]pub unsafe fn bitop_avx512(lhs:&mut[Block],rhs:&[Block]){for i in 0..lhs.len(){let lhs_value=_mm512_load_si512(lhs.as_ptr().add(i).cast());let rhs_value=_mm512_load_si512(rhs.as_ptr().add(i).cast());let value=match OP{BIT_AND=>_mm512_and_si512(lhs_value,rhs_value),BIT_OR=>_mm512_or_si512(lhs_value,rhs_value),BIT_XOR=>_mm512_xor_si512(lhs_value,rhs_value),_=>unreachable!(),};_mm512_store_si512(lhs.as_mut_ptr().add(i).cast(),value);}}#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]pub unsafe fn count_ones_avx2(bits:&[Block])->u64{let table=_mm256_set_epi64x(0x0403_0302_0302_0201,0x0302_0201_0201_0100,0x0403_0302_0302_0201,0x0302_0201_0201_0100);let low_mask=_mm256_set1_epi8(0x0f);let zero=_mm256_setzero_si256();let mut sum=zero;let ptr=bits.as_ptr().cast::<__m256i>();let mut i=0;while i+16<=bits.len()*2{let mut counts=zero;for offset in 0..16{let value=_mm256_load_si256(ptr.add(i+offset));let low=_mm256_shuffle_epi8(table,_mm256_and_si256(value,low_mask));let high=_mm256_shuffle_epi8(table,_mm256_and_si256(_mm256_srli_epi16::<4>(value),low_mask));counts=_mm256_add_epi8(counts,_mm256_add_epi8(low,high));}sum=_mm256_add_epi64(sum,_mm256_sad_epu8(counts,zero));i+=16;}while i(value),low_mask));sum=_mm256_add_epi64(sum,_mm256_sad_epu8(_mm256_add_epi8(low,high),zero));i+=1;}let mut lanes=[0;4];_mm256_storeu_si256(lanes.as_mut_ptr().cast(),sum);lanes.into_iter().sum()}#[target_feature(enable="avx512f,avx512vpopcntdq")]#[cfg(target_arch="x86_64")]pub unsafe fn count_ones_avx512(bits:&[Block])->u64{let mut sum=_mm512_setzero_si512();for i in 0..bits.len(){let value=_mm512_load_si512(bits.as_ptr().add(i).cast());sum=_mm512_add_epi64(sum,_mm512_popcnt_epi64(value));}let mut lanes=[0;8];_mm512_storeu_si512(lanes.as_mut_ptr().cast(),sum);lanes.into_iter().sum()}#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]pub unsafe fn none_avx2(bits:&[Block])->bool{let ptr=bits.as_ptr().cast::<__m256i>();for i in 0..bits.len()*2{let value=_mm256_load_si256(ptr.add(i));if _mm256_testz_si256(value,value)==0{return false;}}true}#[target_feature(enable="avx512f")]#[cfg(target_arch="x86_64")]pub unsafe fn none_avx512(bits:&[Block])->bool{let mut i=0;while i+4<=bits.len(){let value=_mm512_or_si512(_mm512_or_si512(_mm512_load_si512(bits.as_ptr().add(i).cast()),_mm512_load_si512(bits.as_ptr().add(i+1).cast())),_mm512_or_si512(_mm512_load_si512(bits.as_ptr().add(i+2).cast()),_mm512_load_si512(bits.as_ptr().add(i+3).cast())));if _mm512_test_epi64_mask(value,value)!=0{return false;}i+=4;}while ibool{let ones=_mm256_set1_epi64x(-1);let ptr=bits.as_ptr().cast::<__m256i>();for i in 0..bits.len(){let value=_mm256_and_si256(_mm256_load_si256(ptr.add(i*2)),_mm256_load_si256(ptr.add(i*2+1)));if _mm256_movemask_epi8(_mm256_cmpeq_epi64(value,ones))!=-1{return false;}}true}#[target_feature(enable="avx512f")]#[cfg(target_arch="x86_64")]pub unsafe fn all_avx512(bits:&[Block])->bool{let ones=_mm512_set1_epi64(-1);let mut i=0;while i+4<=bits.len(){let value=_mm512_and_si512(_mm512_and_si512(_mm512_load_si512(bits.as_ptr().add(i).cast()),_mm512_load_si512(bits.as_ptr().add(i+1).cast())),_mm512_and_si512(_mm512_load_si512(bits.as_ptr().add(i+2).cast()),_mm512_load_si512(bits.as_ptr().add(i+3).cast())));if _mm512_cmpeq_epi64_mask(value,ones)!=u8::MAX{return false;}i+=4;}while i(bits:&mut[u64],rhs:usize){let word_shift=rhs>>6;let bit_shift=rhs&63;let lower=word_shift+usize::from(bit_shift!=0);let mut end=bits.len();let count=_mm_cvtsi64_si128(bit_shift as i64);while end>=lower+4{let start=end-4;let value=_mm256_loadu_si256(bits.as_ptr().add(start-word_shift).cast());let mut value=if bit_shift==0{value}else{_mm256_or_si256(_mm256_sll_epi64(value,count),_mm256_srl_epi64(_mm256_loadu_si256(bits.as_ptr().add(start-word_shift-1).cast()),_mm_cvtsi64_si128((64-bit_shift)as i64)))};if OR_ASSIGN{value=_mm256_or_si256(value,_mm256_loadu_si256(bits.as_ptr().add(start).cast()));}_mm256_storeu_si256(bits.as_mut_ptr().add(start).cast(),value);end=start;}for i in(lower..end).rev(){let source=i-word_shift;let value=if bit_shift==0{bits[source]}else{(bits[source]<>(64-bit_shift))};if OR_ASSIGN{bits[i]|=value;}else{bits[i]=value;}}if bit_shift!=0{if OR_ASSIGN{bits[word_shift]|=bits[0]<(bits:&mut[u64],rhs:usize){let word_shift=rhs>>6;let bit_shift=rhs&63;let upper=bits.len()-word_shift-usize::from(bit_shift!=0);let mut start=0;let count=_mm_cvtsi64_si128(bit_shift as i64);while start+4<=upper{let value=_mm256_loadu_si256(bits.as_ptr().add(start+word_shift).cast());let mut value=if bit_shift==0{value}else{_mm256_or_si256(_mm256_srl_epi64(value,count),_mm256_sll_epi64(_mm256_loadu_si256(bits.as_ptr().add(start+word_shift+1).cast()),_mm_cvtsi64_si128((64-bit_shift)as i64)))};if OR_ASSIGN{value=_mm256_or_si256(value,_mm256_loadu_si256(bits.as_ptr().add(start).cast()));}_mm256_storeu_si256(bits.as_mut_ptr().add(start).cast(),value);start+=4;}for i in start..upper{let source=i+word_shift;let value=if bit_shift==0{bits[source]}else{(bits[source]>>bit_shift)|(bits[source+1]<<(64-bit_shift))};if OR_ASSIGN{bits[i]|=value;}else{bits[i]=value;}}if bit_shift!=0{if OR_ASSIGN{bits[bits.len()-word_shift-1]|=bits[bits.len()-1]>>bit_shift;}else{bits[bits.len()-word_shift-1]=bits[bits.len()-1]>>bit_shift;}}if!OR_ASSIGN{let end=bits.len()-word_shift;bits[end..].fill(0);}}#[target_feature(enable="avx512f,avx512vbmi2")]#[cfg(target_arch="x86_64")]pub unsafe fn shift_left_avx512(bits:&mut[Block],rhs:usize){let block_shift=rhs>>9;let word_shift=rhs>>6&7;let bit_shift=rhs&63;if word_shift|bit_shift==0{for i in(block_shift..bits.len()).rev(){let mut value=_mm512_load_si512(bits.as_ptr().add(i-block_shift).cast());if OR_ASSIGN{value=_mm512_or_si512(value,_mm512_load_si512(bits.as_ptr().add(i).cast()));}_mm512_store_si512(bits.as_mut_ptr().add(i).cast(),value);}}else{let zero=_mm512_setzero_si512();let indices=_mm512_setr_epi64(0,1,2,3,4,5,6,7);let current_indices=_mm512_add_epi64(indices,_mm512_set1_epi64((8-word_shift)as i64));let previous_indices=_mm512_sub_epi64(current_indices,_mm512_set1_epi64(1));let count=_mm512_set1_epi64(bit_shift as i64);for i in(block_shift..bits.len()).rev(){let source=i-block_shift;let previous=if source==0{zero}else{_mm512_load_si512(bits.as_ptr().add(source-1).cast())};let current=_mm512_load_si512(bits.as_ptr().add(source).cast());let value=if word_shift==0{current}else{_mm512_permutex2var_epi64(previous,current_indices,current)};let mut value=if bit_shift==0{value}else{_mm512_shldv_epi64(value,_mm512_permutex2var_epi64(previous,previous_indices,current),count)};if OR_ASSIGN{value=_mm512_or_si512(value,_mm512_load_si512(bits.as_ptr().add(i).cast()));}_mm512_store_si512(bits.as_mut_ptr().add(i).cast(),value);}}if!OR_ASSIGN{bits[..block_shift].fill(Block::default());}}#[target_feature(enable="avx512f,avx512vbmi2")]#[cfg(target_arch="x86_64")]pub unsafe fn shift_right_avx512(bits:&mut[Block],rhs:usize){let block_shift=rhs>>9;let word_shift=rhs>>6&7;let bit_shift=rhs&63;let remaining=bits.len()-block_shift;if word_shift|bit_shift==0{for i in 0..remaining{let mut value=_mm512_load_si512(bits.as_ptr().add(i+block_shift).cast());if OR_ASSIGN{value=_mm512_or_si512(value,_mm512_load_si512(bits.as_ptr().add(i).cast()));}_mm512_store_si512(bits.as_mut_ptr().add(i).cast(),value);}}else{let zero=_mm512_setzero_si512();let indices=_mm512_setr_epi64(0,1,2,3,4,5,6,7);let current_indices=_mm512_add_epi64(indices,_mm512_set1_epi64(word_shift as i64));let next_indices=_mm512_add_epi64(current_indices,_mm512_set1_epi64(1));let count=_mm512_set1_epi64(bit_shift as i64);for i in 0..remaining{let source=i+block_shift;let current=_mm512_load_si512(bits.as_ptr().add(source).cast());let next=if source+1==bits.len(){zero}else{_mm512_load_si512(bits.as_ptr().add(source+1).cast())};let value=if word_shift==0{current}else{_mm512_permutex2var_epi64(current,current_indices,next)};let mut value=if bit_shift==0{value}else{_mm512_shrdv_epi64(value,_mm512_permutex2var_epi64(current,next_indices,next),count)};if OR_ASSIGN{value=_mm512_or_si512(value,_mm512_load_si512(bits.as_ptr().add(i).cast()));}_mm512_store_si512(bits.as_mut_ptr().add(i).cast(),value);}}if!OR_ASSIGN{bits[remaining..].fill(Block::default());}}}} // 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 Matrixwhere R:SemiRing{pub shape:(usize,usize),pub data:Vec>,_marker:PhantomDataR>}implDebug for Matrixwhere R:SemiRing{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()}}implClone for Matrixwhere R:SemiRing{fn clone(&self)->Self{Self{shape:self.shape,data:self.data.clone(),_marker:self._marker}}}implPartialEq for Matrixwhere R:SemiRing{fn eq(&self,other:&Self)->bool{self.shape==other.shape&&self.data==other.data}}implEq for Matrixwhere R:SemiRing{}implMatrixwhere 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>)->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(&self,mut f:F)->Matrixwhere S:SemiRing,F:FnMut(&R::T)->S::T{Matrix::::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(&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 SystemOfLinearEquationsSolutionwhere R:Field{pub particular:Vec,pub basis:Vec>}implMatrixwhere R:Field{fn eliminate(&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=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(&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::().0}pub fn determinant(&mut self)->R::T{assert_eq!(self.shape.0,self.shape.1);self.eliminate::().1}pub fn solve_system_of_linear_equations(&self,b:&[R::T])->Option>{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::().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>=(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>{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{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>=(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()}}implIndexfor Matrixwhere R:SemiRing{type Output=Vec;fn index(&self,index:usize)->&Self::Output{&self.data[index]}}implIndexMutfor Matrixwhere R:SemiRing{fn index_mut(&mut self,index:usize)->&mut Self::Output{&mut self.data[index]}}implIndex<(usize,usize)>for Matrixwhere R:SemiRing{type Output=R::T;fn index(&self,index:(usize,usize))->&Self::Output{&self.data[index.0][index.1]}}implIndexMut<(usize,usize)>for Matrixwhere 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$imp_assign for Matrixwhere R:SemiRing,$($($clauses)*)?{fn$method_assign(&mut self,rhs:Self){self.pairwise_assign(&rhs,|a,b|R::$method_assign(a,b));}}impl$imp_assign<&Matrix>for Matrixwhere R:SemiRing,$($($clauses)*)?{fn$method_assign(&mut self,rhs:&Self){self.pairwise_assign(rhs,|a,b|R::$method_assign(a,b));}}impl$imp for Matrixwhere R:SemiRing,$($($clauses)*)?{type Output=Matrix;fn$method(mut self,rhs:Self)->Self::Output{self.$method_assign(rhs);self}}impl$imp<&Matrix>for Matrixwhere R:SemiRing,$($($clauses)*)?{type Output=Matrix;fn$method(mut self,rhs:&Self)->Self::Output{self.$method_assign(rhs);self}}impl$imp>for&Matrixwhere R:SemiRing,$($($clauses)*)?{type Output=Matrix;fn$method(self,mut rhs:Matrix)->Self::Output{rhs.pairwise_assign(self,|a,b|*a=R::$method(b,a));rhs}}impl$imp<&Matrix>for&Matrixwhere R:SemiRing,$($($clauses)*)?{type Output=Matrix;fn$method(self,rhs:&Matrix)->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]);implMul for Matrixwhere R:SemiRing{type Output=Matrix;fn mul(self,rhs:Self)->Self::Output{(&self).mul(&rhs)}}implMul<&Matrix>for Matrixwhere R:SemiRing{type Output=Matrix;fn mul(self,rhs:&Matrix)->Self::Output{(&self).mul(rhs)}}implMul>for&Matrixwhere R:SemiRing{type Output=Matrix;fn mul(self,rhs:Matrix)->Self::Output{self.mul(&rhs)}}implMul<&Matrix>for&Matrixwhere R:SemiRing{type Output=Matrix;fn mul(self,rhs:&Matrix)->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(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(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(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::(&a[a11..],&a[a22..],s1,k,stride_a,stride_a);add_block::(&b[b11..],&b[b22..],s2,w,stride_b,stride_b);strassen_rec::(s1,s2,m1,(h,k,w),k,w);add_block::(&a[a21..],&a[a22..],s1,k,stride_a,stride_a);strassen_rec::(s1,&b[b11..],m2,(h,k,w),k,stride_b);sub_block::(&b[b12..],&b[b22..],s2,w,stride_b,stride_b);strassen_rec::(&a[a11..],s2,m3,(h,k,w),stride_a,w);sub_block::(&b[b21..],&b[b11..],s2,w,stride_b,stride_b);strassen_rec::(&a[a22..],s2,m4,(h,k,w),stride_a,w);add_block::(&a[a11..],&a[a12..],s1,k,stride_a,stride_a);strassen_rec::(s1,&b[b22..],m5,(h,k,w),k,stride_b);sub_block::(&a[a21..],&a[a11..],s1,k,stride_a,stride_a);add_block::(&b[b11..],&b[b12..],s2,w,stride_b,stride_b);strassen_rec::(s1,s2,m6,(h,k,w),k,w);sub_block::(&a[a12..],&a[a22..],s1,k,stride_a,stride_a);add_block::(&b[b21..],&b[b22..],s2,w,stride_b,stride_b);strassen_rec::(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);}}implMatrixwhere R:Ring{pub fn mul_strassen(&self,rhs:&Matrix)->Matrix{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::(&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}}implMulAssign<&R::T>for Matrixwhere 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);}}}}implNeg for Matrixwhere R:SemiRing{type Output=Self;fn neg(self)->Self::Output{self.map(|x|R::neg(x))}}implNeg for&Matrixwhere R:SemiRing{type Output=Matrix;fn neg(self)->Self::Output{self.map(|x|R::neg(x))}}implMatrixwhere 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}}implMatrixwhere 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 MultiplicativeOperationwhere T:Clone+One+Mul{_marker:PhantomDataT>}implMagma for MultiplicativeOperationwhere T:Clone+One+Mul{type T=T;#[inline]fn operate(x:&Self::T,y:&Self::T)->Self::T{x.clone()*y.clone()}}implUnital for MultiplicativeOperationwhere T:Clone+One+Mul{#[inline]fn unit()->Self::T{One::one()}}implAssociative for MultiplicativeOperationwhere T:Clone+One+Mul{}implCommutative for MultiplicativeOperationwhere T:Clone+One+Mul{}implInvertible for MultiplicativeOperationwhere T:Clone+One+Mul+Div{#[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: Xorshift #[cfg_attr(any(),rustfmt::skip)]pub use self::xorshift::Xorshift;#[cfg_attr(any(),rustfmt::skip)]mod xorshift{use std::{collections::hash_map::RandomState,hash::{BuildHasher,Hasher},iter::repeat_with};#[derive(Clone,Debug)]pub struct Xorshift{y:u64}impl Xorshift{pub fn new_with_seed(seed:u64)->Self{Xorshift{y:seed}}pub fn new()->Self{Xorshift::new_with_seed(RandomState::new().build_hasher().finish())}pub fn rand64(&mut self)->u64{self.y^=self.y<<5;self.y^=self.y>>17;self.y^=self.y<<11;self.y}pub fn rand(&mut self,k:u64)->u64{self.rand64()%k}pub fn rands(&mut self,k:u64,n:usize)->Vec{repeat_with(| |self.rand(k)).take(n).collect()}pub fn randf(&mut self)->f64{const UPPER_MASK:u64=0x3FF0_0000_0000_0000;const LOWER_MASK:u64=0x000F_FFFF_FFFF_FFFF;let x=self.rand64();let tmp=UPPER_MASK|(x&LOWER_MASK);let result:f64=f64::from_bits(tmp);f64::from_bits(f64::to_bits(result-1.0)^(x>>63))}pub fn gen_bool(&mut self,p:f64)->bool{self.randf()(&mut self,slice:&mut[T]){let mut n=slice.len();while n>1{let i=self.rand(n as _)as usize;n-=1;slice.swap(i,n);}}}impl Default for Xorshift{fn default()->Self{Xorshift::new_with_seed(0x2b99_2ddf_a232_49d6)}}} // codesnip-guard: _codesnip_and_4d6174726978_636f64696e67 #[cfg_attr(any(),rustfmt::skip)]implSerdeByteStr for Matrixwhere R:SemiRing{fn serialize(&self,buf:&mut Vec){self.data.serialize(buf);}fn deserialize(iter:&mut I)->Self where I:Iterator{Self::from_vec(Vec::deserialize(iter))}} // 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{}implSemiGroup 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;fn bits(self)->Self::Iter;}pub trait SignedExpBits{type T:ExpBits;fn neg_and_bits(self)->(bool,Self::T);}pub struct Bits{n:T}macro_rules!impl_exp_bits_for_uint{($($t:ty)*)=>{$(impl Iterator for Bits<$t>{type Item=bool;fn next(&mut self)->Option{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(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(iter:I)->Self::T where I:IntoIterator{let mut iter=iter.into_iter();if let Some(item)=iter.next(){iter.fold(item,|acc,x|Self::operate(&acc,&x))}else{Self::unit()}}}implMonoid 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(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}}}implGroup 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{}implAbelianMonoid for M where M:Monoid+Commutative{}#[doc=" commutative group"]pub trait AbelianGroup:Group+Commutative{}implAbelianGroup 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{}implIdempotentMonoid 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: avx_helper #[cfg_attr(any(),rustfmt::skip)]pub use self::avx_helper::{SimdBackend,avx512_enabled,avx512_supported,disable_avx512,enable_avx512,simd_backend};#[cfg_attr(any(),rustfmt::skip)]mod avx_helper{use std::sync::atomic::{AtomicBool,Ordering};#[derive(Copy,Clone,Debug,Eq,PartialEq)]pub enum SimdBackend{Scalar,Avx2,Avx512}static AVX512_ENABLED:AtomicBool=AtomicBool::new(true);#[inline]pub fn disable_avx512(){AVX512_ENABLED.store(false,Ordering::Relaxed);}#[inline]pub fn enable_avx512(){AVX512_ENABLED.store(true,Ordering::Relaxed);}#[inline]pub fn avx512_enabled()->bool{AVX512_ENABLED.load(Ordering::Relaxed)}#[inline]pub fn avx512_supported()->bool{#[cfg(any(target_arch="x86",target_arch="x86_64"))]return is_x86_feature_detected!("avx512f")&&is_x86_feature_detected!("avx512dq")&&is_x86_feature_detected!("avx512cd")&&is_x86_feature_detected!("avx512bw")&&is_x86_feature_detected!("avx512vl");#[cfg(not(any(target_arch="x86",target_arch="x86_64")))]false}#[inline]pub fn simd_backend()->SimdBackend{#[cfg(any(target_arch="x86",target_arch="x86_64"))]{if avx512_enabled()&&avx512_supported(){return SimdBackend::Avx512;}if is_x86_feature_detected!("avx2"){return SimdBackend::Avx2;}}SimdBackend::Scalar}#[macro_export]macro_rules!avx_helper{(@dispatch$backend:path,$kind:ident;$avx512:expr,$avx2:expr,$scalar:expr)=>{{#[cfg(target_arch="x86_64")]{match$backend(){$kind::Avx512=>$avx512,$kind::Avx2=>$avx2,$kind::Scalar=>$scalar,}}#[cfg(not(target_arch="x86_64"))]$scalar}};(@dispatch_avx2_fma$avx2:expr,$scalar:expr)=>{{#[cfg(target_arch="x86_64")]{if is_x86_feature_detected!("avx2")&&is_x86_feature_detected!("fma"){$avx2}else{$scalar}}#[cfg(not(target_arch="x86_64"))]$scalar}};(@avx512$(#[$meta:meta])*$vis:vis fn$name:ident$(<$($T:ident),+>)?($($i:ident:$t:ty),*)->$ret:ty where[$($clauses:tt)*]$body:block)=>{$(#[$meta])*$vis fn$name$(<$($T)*>)?($($i:$t),*)->$ret where$($clauses)*{if$crate::avx512_supported(){$crate::avx_helper!(@def_avx512 fn avx512$(<$($T)*>)?($($i:$t),*)->$ret where[$($clauses)*]$body);unsafe{avx512$(::<$($T),*>)?($($i),*)}}else if is_x86_feature_detected!("avx2"){$crate::avx_helper!(@def_avx2 fn avx2$(<$($T)*>)?($($i:$t),*)->$ret where[$($clauses)*]$body);unsafe{avx2$(::<$($T),*>)?($($i),*)}}else{$body}}};(@avx2$(#[$meta:meta])*$vis:vis fn$name:ident$(<$($T:ident),+>)?($($i:ident:$t:ty),*)->$ret:ty where[$($clauses:tt)*]$body:block)=>{$(#[$meta])*$vis fn$name$(<$($T)*>)?($($i:$t),*)->$ret where$($clauses)*{if is_x86_feature_detected!("avx2"){$crate::avx_helper!(@def_avx2 fn avx2$(<$($T)*>)?($($i:$t),*)->$ret where[$($clauses)*]$body);unsafe{avx2$(::<$($T),*>)?($($i),*)}}else{$body}}};(@def_avx512 fn$name:ident$(<$($T:ident),+>)?($($args:tt)*)->$ret:ty where[$($clauses:tt)*]$body:block)=>{#[target_feature(enable="avx512f,avx512dq,avx512cd,avx512bw,avx512vl")]unsafe fn$name$(<$($T)*>)?($($args)*)->$ret where$($clauses)*$body};(@def_avx2 fn$name:ident$(<$($T:ident),+>)?($($args:tt)*)->$ret:ty where[$($clauses:tt)*]$body:block)=>{#[target_feature(enable="avx2")]unsafe fn$name$(<$($T)*>)?($($args)*)->$ret where$($clauses)*$body};(@$tag:ident$(#[$meta:meta])*$vis:vis fn$name:ident$(<$($T:ident),+>)?($($args:tt)*)->$ret:ty$body:block)=>{$crate::avx_helper!(@$tag$(#[$meta])*$vis fn$name$(<$($T)*>)?($($args)*)->$ret where[]$body);};(@$tag:ident$(#[$meta:meta])*$vis:vis fn$name:ident$(<$($T:ident),+>)?($($args:tt)*)$($t:tt)*)=>{$crate::avx_helper!(@$tag$(#[$meta])*$vis fn$name$(<$($T)*>)?($($args)*)->()$($t)*);};($($t:tt)*)=>{::std::compile_error!($($t)*);}}} // codesnip-guard: coding #[cfg_attr(any(),rustfmt::skip)]pub use self::coding::{SerdeByteStr,unescape};#[cfg_attr(any(),rustfmt::skip)]mod coding{use std::{cmp::Reverse,collections::{BTreeMap,BTreeSet,BinaryHeap,HashMap,HashSet,VecDeque},hash::Hash,iter::repeat_with,mem::size_of};pub fn unescape(bytes:&[u8])->Vec{let mut buf=Vec::with_capacity(bytes.len());let mut escape=false;for b in bytes{match(escape,*b){(true,b'n')=>{buf.push(b'\n');escape=false;}(true,b'r')=>{buf.push(b'\r');escape=false;}(true,b't')=>{buf.push(b'\t');escape=false;}(true,b'\\')=>{buf.push(b'\\');escape=false;}(true,b'0')=>{buf.push(b'\0');escape=false;}(true,b'"')=>{buf.push(b'"');escape=false;}(false,b'\\')=>{escape=true;}(_,b)=>{buf.push(b);escape=false;}}}buf}fn to_bytestring(bytes:&[u8])->String{let mut count=[0usize;128];let mut c=0u8;let mut w=0u32;for&b in bytes{w+=1;count[(c|(b>>w))as usize]+=1;if w==7{count[(b&0b1111111u8)as usize]+=1;c=0;w=0;}else{c=(b<<(7-w))&0b1111111u8;}}if w>0{count[c as usize]+=1;c=0;w=0;}let mut salt=0u8;let mut extra_min=!0usize;for s in 0u8..128{let mut extra=0usize;for&b in b"\n\r\t\\\0\""{extra+=count[(b^s)as usize];extra+=(b==s)as usize;}if extra()+1;let mut buf=String::with_capacity(cap);macro_rules!escape_branch{($e:literal)=>{{buf.push('\\');buf.push($e);}};}macro_rules!escape{($b:expr)=>{match($b)^salt{b'\n'=>escape_branch!('n'),b'\r'=>escape_branch!('r'),b'\t'=>escape_branch!('t'),b'\\'=>escape_branch!('\\'),b'\0'=>escape_branch!('0'),b'"'=>escape_branch!('"'),b=>buf.push(b as char),}};}escape!(0);for&b in bytes{w+=1;escape!(c|(b>>w));if w==7{escape!(b&0b1111111u8);c=0;w=0;}else{c=(b<<(7-w))&0b1111111u8;}}if w>0{escape!(c);}assert_eq!(cap,buf.len());buf}fn from_bytestring(bytes:&[u8])->Vec{assert!(!bytes.is_empty());let cap=(bytes.len()-1)*7/8;let mut buf=Vec::with_capacity(cap);let salt=bytes[0];let bytes=&bytes[1..];let mut c=0u8;let mut w=0u32;for&b in bytes{let b=b^salt;if w==0{c=b<<1;w=7;}else{w-=1;buf.push(c|(b>>w));c=if w>0{b<<(8-w)}else{0};}}assert_eq!(cap,buf.len());buf}#[derive(Debug,PartialEq,Eq,PartialOrd,Ord)]enum HuffmanTree{Leaf(u8),Node(Box,Box)}#[derive(Debug)]struct BitWriter{bytes:Vec,last:u8,w:u32}impl Default for BitWriter{fn default()->Self{Self{bytes:Default::default(),last:Default::default(),w:8}}}impl BitWriter{fn push_bit(&mut self,b:bool){self.w-=1;self.last|=(b as u8)<>(8-self.w)));if self.w<8{self.last=b<>c<>c)as u8);}self.last=0;self.w=8-t;if t>0{self.last=(b as u8)<Vec{if self.w<8{self.bytes.push(self.last);}self.bytes}}#[derive(Debug)]struct BitReader<'a>{bytes:&'a[u8],pos:u32}impl<'a>BitReader<'a>{fn new(bytes:&'a[u8])->Self{Self{bytes,pos:0}}fn read_bit(&mut self)->bool{let b=(self.bytes[0]>>(7-self.pos))&1==1;self.pos+=1;if self.pos==8{self.pos=0;self.bytes=&self.bytes[1..];}b}fn read_u8(&mut self)->u8{let b=self.bytes[0]<>(8-self.pos))}}}fn huffman_coding(bytes:&[u8])->Vec{fn make_table(t:&HuffmanTree,code:u64,len:u32,table:&mut[(u64,u32)]){match t{HuffmanTree::Leaf(i)=>{table[*i as usize]=(code,len.max(1));}HuffmanTree::Node(l,r)=>{make_table(l,code<<1,len+1,table);make_table(r,(code<<1)|1,len+1,table);}}}fn output_tree(t:&HuffmanTree,buf:&mut BitWriter){match t{HuffmanTree::Leaf(i)=>{buf.push_bit(false);buf.push_u8(*i);}HuffmanTree::Node(l,r)=>{buf.push_bit(true);output_tree(l,buf);output_tree(r,buf);}}}let mut freq=[0usize;256];for&b in&bytes.len().to_le_bytes(){freq[b as usize]+=1;}for&b in bytes{freq[b as usize]+=1;}let mut heap=BinaryHeap::new();for(i,&f)in freq.iter().enumerate(){if f>0{heap.push(Reverse((f,0usize,HuffmanTree::Leaf(i as _))));}}let t=if heap.is_empty(){HuffmanTree::Node(Box::new(HuffmanTree::Leaf(0)),Box::new(HuffmanTree::Leaf(0)))}else{loop{let Reverse((f,c,t))=heap.pop().unwrap();if let Some(Reverse((ff,cc,tt)))=heap.pop(){heap.push(Reverse((f+ff,c.max(cc)+1,HuffmanTree::Node(Box::new(t),Box::new(tt)))));}else{break t;}}};let mut table=vec![(0u64,0u32);256];make_table(&t,0,0,&mut table);let mut buf=BitWriter::default();output_tree(&t,&mut buf);for&b in&bytes.len().to_le_bytes(){let(x,y)=table[b as usize];buf.push_u64(x,y);}for&b in bytes{let(x,y)=table[b as usize];buf.push_u64(x,y);}buf.into_inner()}fn huffman_decoding(bytes:&[u8])->Vec{fn read_tree(reader:&mut BitReader)->HuffmanTree{if reader.read_bit(){HuffmanTree::Node(Box::new(read_tree(reader)),Box::new(read_tree(reader)))}else{HuffmanTree::Leaf(reader.read_u8())}}fn decode(mut t:&HuffmanTree,reader:&mut BitReader)->u8{loop{match t{HuffmanTree::Leaf(i)=>break*i,HuffmanTree::Node(l,r)=>t=if reader.read_bit(){r}else{l},}}}let mut reader=BitReader::new(bytes);let t=read_tree(&mut reader);const C:usize=size_of::();let mut size:[u8;C]=[0u8;C];for b in&mut size{*b=decode(&t,&mut reader);}let size=usize::from_le_bytes(size);let mut buf=vec![];for i in 0..size{buf.push(decode(&t,&mut reader));if i<10{}}buf}pub trait SerdeByteStr{fn serialize(&self,buf:&mut Vec);fn deserialize(iter:&mut I)->Self where I:Iterator;fn serialize_bytestr(&self)->String{let mut bytes=vec![];self.serialize(&mut bytes);let bytes=huffman_coding(&bytes);to_bytestring(&bytes)}fn deserialize_from_bytes(bytes:&[u8])->Self where Self:Sized{let bytes=from_bytestring(bytes);let bytes=huffman_decoding(&bytes);Self::deserialize(&mut bytes.iter().cloned())}}impl SerdeByteStr for bool{fn serialize(&self,buf:&mut Vec){(*self as u8).serialize(buf)}fn deserialize(iter:&mut I)->Self where I:Iterator{iter.next().unwrap()!=0}}macro_rules!impl_serdebytestr_num{($($t:ty)*)=>{$(impl SerdeByteStr for$t{fn serialize(&self,buf:&mut Vec){buf.extend(self.to_le_bytes().iter());}fn deserialize(iter:&mut I)->Self where I:Iterator,{const C:usize=size_of::<$t>();let mut bytes:[u8;C]=[0u8;C];for(b,i)in bytes.iter_mut().zip(iter){*b=i;}<$t>::from_le_bytes(bytes)}})*};}impl_serdebytestr_num!(u8 u16 u32 u64 u128 usize i8 i16 i32 i64 i128 isize f32 f64);impl SerdeByteStr for char{fn serialize(&self,buf:&mut Vec){(*self as u32).serialize(buf)}fn deserialize(iter:&mut I)->Self where I:Iterator{char::from_u32(u32::deserialize(iter)).expect("invalid Unicode scalar value")}}impl SerdeByteStr for(){fn serialize(&self,_buf:&mut Vec){}fn deserialize(_iter:&mut I)->Self where I:Iterator{}}macro_rules!impl_serdebytestr_tuple{(@impl$($A:ident$a:ident)*)=>{impl<$($A,)*>SerdeByteStr for($($A,)*)where$($A:SerdeByteStr),*{fn serialize(&self,buf:&mut Vec){let($($a,)*)=self;$(SerdeByteStr::serialize($a,buf));*}fn deserialize(iter:&mut Iter)->Self where Iter:Iterator,{($(<$A as SerdeByteStr>::deserialize(iter),)*)}}};(@inc,$B:ident$b:ident$($C:ident$c:ident)*)=>{impl_serdebytestr_tuple!(@inc$B$b,$($C$c)*);};(@inc$($A:ident$a:ident)*,$B:ident$b:ident$($C:ident$c:ident)*)=>{impl_serdebytestr_tuple!(@impl$($A$a)*);impl_serdebytestr_tuple!(@inc$($A$a)*$B$b,$($C$c)*);};(@inc$($A:ident$a:ident)*,)=>{impl_serdebytestr_tuple!(@impl$($A$a)*);};($($t:tt)*)=>{impl_serdebytestr_tuple!(@inc,$($t)*);};}impl_serdebytestr_tuple!(A a B b C c D d E e F f G g H h I i J j K k);implSerdeByteStr for Optionwhere T:SerdeByteStr{fn serialize(&self,buf:&mut Vec){self.is_some().serialize(buf);if let Some(x)=self{x.serialize(buf);}}fn deserialize(iter:&mut I)->Self where I:Iterator{if bool::deserialize(iter){Some(T::deserialize(iter))}else{None}}}impl SerdeByteStr for String{fn serialize(&self,buf:&mut Vec){let bytes=self.bytes();bytes.len().serialize(buf);for x in bytes{x.serialize(buf);}}fn deserialize(iter:&mut I)->Self where I:Iterator{let n=usize::deserialize(iter);let bytes=match iter.size_hint(){(lower,Some(upper))if lower==upper&&lower>=n=>iter.by_ref().take(n).collect(),_=>repeat_with(| |u8::deserialize(iter)).take(n).collect(),};String::from_utf8(bytes).expect("invalid UTF-8")}}macro_rules!impl_serdebytestr_seq{($([$($g:ident)*]$t:ty$(where[$($tt:tt)*])?),*$(,)?)=>{$(impl<$($g),*>SerdeByteStr for$t where$($g:SerdeByteStr,)*$($($tt)*)?{fn serialize(&self,buf:&mut Vec){self.len().serialize(buf);for x in self{x.serialize(buf);}}fn deserialize(iter:&mut I)->Self where I:Iterator,{let n=usize::deserialize(iter);repeat_with(||SerdeByteStr::deserialize(iter)).take(n).collect()}})*};(@kv$([$($g:ident)*]$t:ty$(where[$($tt:tt)*])?),*$(,)?)=>{$(impl<$($g),*>SerdeByteStr for$t where$($g:SerdeByteStr,)*$($($tt)*)?{fn serialize(&self,buf:&mut Vec){self.len().serialize(buf);for(k,v)in self{k.serialize(buf);v.serialize(buf);}}fn deserialize(iter:&mut I)->Self where I:Iterator,{let n=usize::deserialize(iter);repeat_with(||SerdeByteStr::deserialize(iter)).take(n).collect()}})*};}impl_serdebytestr_seq!([T]Vec,[T]VecDeque,[T]BinaryHeapwhere[T:Ord],[T]BTreeSetwhere[T:Ord],[T]HashSetwhere[T:Eq+Hash],);impl_serdebytestr_seq!(@kv[K V]BTreeMapwhere[K:Ord],[K V]HashMapwhere[K:Eq+Hash],);} // codesnip-guard: random_generator #[cfg_attr(any(),rustfmt::skip)]pub use self::random_generator::{NotEmptySegment,RandIter,RandRange,RandomSpec,Spec,WeightedSampler,WithEmptySegment,from_fn,spec};#[cfg_attr(any(),rustfmt::skip)]mod random_generator{use super::Xorshift;use std::{marker::PhantomData,mem::swap,ops::{Bound,Range,RangeFrom,RangeFull,RangeInclusive,RangeTo,RangeToInclusive}};#[doc=" Trait for spec of generating random value."]pub trait RandomSpec:Sized{#[doc=" Return a random value."]fn rand(&self,rng:&mut Xorshift)->T;#[doc=" Return an iterator that generates random values."]fn rand_iter(self,rng:&mut Xorshift)->RandIter<'_,T,Self>{RandIter{spec:self,rng,_marker:PhantomData}}}pub struct Spec{generate:F,_marker:PhantomDataT>}#[doc=" Wraps a random specification for composition."]pub fn spec(source:R)->SpecT>where R:RandomSpec{from_fn(move|rng|source.rand(rng))}#[doc=" Calls `generate` with the RNG whenever a value is requested."]pub fn from_fn(generate:F)->Specwhere F:Fn(&mut Xorshift)->T{Spec{generate,_marker:PhantomData}}implSpecwhere F:Fn(&mut Xorshift)->T{#[doc=" Samples the length once, then generates that many elements."]pub fn vec(self,length:R)->Spec,impl Fn(&mut Xorshift)->Vec>where R:RandomSpec{from_fn(move|rng|{let n=rng.random(&length);rng.random_iter(&self).take(n).collect()})}}implRandomSpecfor Specwhere F:Fn(&mut Xorshift)->T{fn rand(&self,rng:&mut Xorshift)->T{(self.generate)(rng)}}impl Xorshift{pub fn random(&mut self,spec:R)->T where R:RandomSpec{spec.rand(self)}pub fn random_iter(&mut self,spec:R)->RandIter<'_,T,R>where R:RandomSpec{spec.rand_iter(self)}}#[derive(Debug)]pub struct RandIter<'r,T,R>where R:RandomSpec{spec:R,rng:&'r mut Xorshift,_marker:PhantomDataT>}implIterator for RandIter<'_,T,R>where R:RandomSpec{type Item=T;fn next(&mut self)->Option{Some(self.spec.rand(self.rng))}fn size_hint(&self)->(usize,Option){(usize::MAX,None)}}macro_rules!impl_random_spec_range_full{($($t:ty)*)=>{$(impl RandomSpec<$t>for RangeFull{fn rand(&self,rng:&mut Xorshift)->$t{rng.rand64()as _}})*};}impl_random_spec_range_full!(u8 u16 u32 u64 usize i8 i16 i32 i64 isize);impl RandomSpecfor RangeFull{fn rand(&self,rng:&mut Xorshift)->u128{((rng.rand64()as u128)<<64)|rng.rand64()as u128}}impl RandomSpecfor RangeFull{fn rand(&self,rng:&mut Xorshift)->i128{rng.random::(..)as i128}}macro_rules!impl_random_spec_ranges{($($u:ident$i:ident)*)=>{$(impl RandomSpec<$u>for Range<$u>{fn rand(&self,rng:&mut Xorshift)->$u{assert!(self.start(..)%len)}}impl RandomSpec<$i>for Range<$i>{fn rand(&self,rng:&mut Xorshift)->$i{assert!(self.start(..)%len)}}impl RandomSpec<$u>for RangeFrom<$u>{fn rand(&self,rng:&mut Xorshift)->$u{let len=($u::MAX-self.start).wrapping_add(1);let x=rng.random::<$u,_>(..);self.start+if len!=0{x%len}else{x}}}impl RandomSpec<$i>for RangeFrom<$i>{fn rand(&self,rng:&mut Xorshift)->$i{let len=($i::MAX.abs_diff(self.start)).wrapping_add(1);let x=rng.random::<$u,_>(..);self.start.wrapping_add_unsigned(if len!=0{x%len}else{x})}}impl RandomSpec<$u>for RangeInclusive<$u>{fn rand(&self,rng:&mut Xorshift)->$u{assert!(self.start()<=self.end());let len=(self.end()-self.start()).wrapping_add(1);let x=rng.random::<$u,_>(..);self.start()+if len!=0{x%len}else{x}}}impl RandomSpec<$i>for RangeInclusive<$i>{fn rand(&self,rng:&mut Xorshift)->$i{assert!(self.start()<=self.end());let len=(self.end().abs_diff(*self.start())).wrapping_add(1);let x=rng.random::<$u,_>(..);self.start().wrapping_add_unsigned(if len!=0{x%len}else{x})}}impl RandomSpec<$u>for RangeTo<$u>{fn rand(&self,rng:&mut Xorshift)->$u{let len=self.end;rng.random::<$u,_>(..)%len}}impl RandomSpec<$i>for RangeTo<$i>{fn rand(&self,rng:&mut Xorshift)->$i{let len=self.end.abs_diff($i::MIN);$i::MIN.wrapping_add_unsigned(rng.random::<$u,_>(..)%len)}}impl RandomSpec<$u>for RangeToInclusive<$u>{fn rand(&self,rng:&mut Xorshift)->$u{let len=(self.end).wrapping_add(1);let x=rng.random::<$u,_>(..);if len!=0{x%len}else{x}}}impl RandomSpec<$i>for RangeToInclusive<$i>{fn rand(&self,rng:&mut Xorshift)->$i{let len=(self.end.abs_diff($i::MIN)).wrapping_add(1);let x=rng.random::<$u,_>(..);$i::MIN.wrapping_add_unsigned(if len!=0{x%len}else{x})}})*};}impl_random_spec_ranges!(u8 i8 u16 i16 u32 i32 u64 i64 u128 i128 usize isize);macro_rules!impl_random_spec_tuple{($($T:ident)*,$($R:ident)*,$($v:ident)*)=>{impl<$($T),*,$($R),*>RandomSpec<($($T,)*)>for($($R,)*)where$($R:RandomSpec<$T>),*{fn rand(&self,rng:&mut Xorshift)->($($T,)*){let($($v,)*)=self;($(($v).rand(rng),)*)}}};}impl_random_spec_tuple!(A,RA,a);impl_random_spec_tuple!(A B,RA RB,a b);impl_random_spec_tuple!(A B C,RA RB RC,a b c);impl_random_spec_tuple!(A B C D,RA RB RC RD,a b c d);impl_random_spec_tuple!(A B C D E,RA RB RC RD RE,a b c d e);impl_random_spec_tuple!(A B C D E F,RA RB RC RD RE RF,a b c d e f);impl_random_spec_tuple!(A B C D E F G,RA RB RC RD RE RF RG,a b c d e f g);impl_random_spec_tuple!(A B C D E F G H,RA RB RC RD RE RF RG RH,a b c d e f g h);impl_random_spec_tuple!(A B C D E F G H I,RA RB RC RD RE RF RG RH RI,a b c d e f g h i);impl_random_spec_tuple!(A B C D E F G H I J,RA RB RC RD RE RF RG RH RI RJ,a b c d e f g h i j);macro_rules!impl_random_spec_primitive{($($t:ty)*)=>{$(impl RandomSpec<$t>for$t{fn rand(&self,_rng:&mut Xorshift)->$t{*self}})*};}impl_random_spec_primitive!(()u8 u16 u32 u64 u128 usize i8 i16 i32 i64 i128 isize bool char);implRandomSpecfor&R where R:RandomSpec{fn rand(&self,rng:&mut Xorshift)->T{>::rand(self,rng)}}implRandomSpecfor&mut R where R:RandomSpec{fn rand(&self,rng:&mut Xorshift)->T{>::rand(self,rng)}}#[derive(Debug,Clone,Copy,PartialEq,Eq,PartialOrd,Ord,Hash)]#[doc=" Left-close Right-open No Empty Segment"]pub struct NotEmptySegment(pub T);implRandomSpec<(usize,usize)>for NotEmptySegmentwhere T:RandomSpec{fn rand(&self,rng:&mut Xorshift)->(usize,usize){let n=rng.random(&self.0)as u64;let k=randint_uniform(rng,n);let l=randint_uniform(rng,n-k)as usize;(l,l+k as usize+1)}}#[derive(Debug,Clone,Copy,PartialEq,Eq,PartialOrd,Ord,Hash)]#[doc=" Left-close Right-open Segment"]pub struct WithEmptySegment(pub T);implRandomSpec<(usize,usize)>for WithEmptySegmentwhere T:RandomSpec{fn rand(&self,rng:&mut Xorshift)->(usize,usize){let n=rng.random(&self.0)as u64;let k=randint_uniform(rng,n+1);let l=randint_uniform(rng,n-k+1)as usize;(l,l+k as usize)}}#[derive(Debug,Clone,Copy,PartialEq,Eq,PartialOrd,Ord,Hash)]pub struct RandRange{data:Q,_marker:PhantomDataT>}implRandRange{pub fn new(data:Q)->Self{Self{data,_marker:PhantomData}}}implRandomSpec<(Bound,Bound)>for RandRangewhere Q:RandomSpec,T:Ord{fn rand(&self,rng:&mut Xorshift)->(Bound,Bound){let mut l=rng.random(&self.data);let mut r=rng.random(&self.data);if l>r{swap(&mut l,&mut r);}(match rng.rand(3){0=>Bound::Excluded(l),1=>Bound::Included(l),_=>Bound::Unbounded,},match rng.rand(3){0=>Bound::Excluded(r),1=>Bound::Included(r),_=>Bound::Unbounded,})}}#[inline]fn randint_uniform(rng:&mut Xorshift,k:u64)->u64{let mut v=rng.rand64();if k>0{v%=k;}v}pub struct WeightedSampler{n:usize,prob:Vec,alias:Vec}impl WeightedSampler{pub fn new(weights:impl IntoIterator)->Self{let mut weights:Vec<_>=weights.into_iter().collect();let n=weights.len();assert!(n>0,"weights must be non-empty");let mut prob=vec![0.0;n];let mut alias=vec![0;n];let mut small=vec![];let mut large=vec![];let sum:f64=weights.iter().sum();assert!(sum>0.0,"sum of weights must be positive");let mut scale=n as f64/sum;if!sum.is_finite()||!scale.is_finite(){let max=weights.iter().fold(0.0f64,|max,&weight|{assert!(weight.is_finite()&&weight>=0.0,"weights must be finite and non-negative");max.max(weight)});let sum:f64=weights.iter_mut().map(|weight|{*weight/=max;*weight}).sum();scale=n as f64/sum;}for(i,weight)in weights.iter_mut().enumerate(){assert!(*weight>=0.0,"weights must be non-negative");*weight*=scale;if*weight<1.0{small.push(i);}else{large.push(i);}}while let Some(&g)=large.last(){let Some(l)=small.pop()else{break;};prob[l]=weights[l];alias[l]=g;weights[g]-=1.0-weights[l];if weights[g]<1.0{large.pop();small.push(g);}}for g in small.into_iter().chain(large).rev(){prob[g]=1.0;alias[g]=g;}Self{n,prob,alias}}}impl RandomSpecfor WeightedSampler{fn rand(&self,rng:&mut Xorshift)->usize{let i=rng.rand(self.n as u64)as usize;if rng.randf(){::std::iter::repeat_with(||$crate::rand_value!(@inner$rng,[]$($t)*))$(.take($len).collect::<::std::vec::Vec<_>>())?};(@array$rng:expr,[$($t:tt)*]$len:expr)=>{{let values:[_;$len]=::std::array::from_fn(|_|$crate::rand_value!(@inner$rng,[]$($t)*));values}};(@group$rng:expr,[])=>{()};(@group$rng:expr,[$($args:tt)*],$($t:tt)*)=>{$crate::rand_value!(@tuple$rng,[$($args)*]$($t)*)};(@tuple$rng:expr,[$([$($args:tt)*])*])=>{($($($args)*,)*)};(@$tag:ident$rng:expr,[[$($args:tt)*]])=>{$($args)*};(@$tag:ident$rng:expr,[$($args:tt)*]($($tuple:tt)*)$(,$($t:tt)*)?)=>{$crate::rand_value!(@$tag$rng,[$($args)*[$crate::rand_value!(@group$rng,[]$($tuple)*)]]$(,$($t)*)?)};(@$tag:ident$rng:expr,[$($args:tt)*][[$($tt:tt)*];const$len:expr]$(,$($t:tt)*)?)=>{$crate::rand_value!(@$tag$rng,[$($args)*[$crate::rand_value!(@array$rng,[[$($tt)*]]$len)]]$(,$($t)*)?)};(@$tag:ident$rng:expr,[$($args:tt)*][[$($tt:tt)*];$len:expr]$(,$($t:tt)*)?)=>{$crate::rand_value!(@$tag$rng,[$($args)*[$crate::rand_value!(@repeat$rng,[[$($tt)*]]$len)]]$(,$($t)*)?)};(@$tag:ident$rng:expr,[$($args:tt)*][($($tt:tt)*);const$len:expr]$(,$($t:tt)*)?)=>{$crate::rand_value!(@$tag$rng,[$($args)*[$crate::rand_value!(@array$rng,[($($tt)*)]$len)]]$(,$($t)*)?)};(@$tag:ident$rng:expr,[$($args:tt)*][($($tt:tt)*);$len:expr]$(,$($t:tt)*)?)=>{$crate::rand_value!(@$tag$rng,[$($args)*[$crate::rand_value!(@repeat$rng,[($($tt)*)]$len)]]$(,$($t)*)?)};(@$tag:ident$rng:expr,[$($args:tt)*][$ty:expr;const$len:expr]$(,$($t:tt)*)?)=>{$crate::rand_value!(@$tag$rng,[$($args)*[$crate::rand_value!(@array$rng,[$ty]$len)]]$(,$($t)*)?)};(@$tag:ident$rng:expr,[$($args:tt)*][$ty:expr;$len:expr]$(,$($t:tt)*)?)=>{$crate::rand_value!(@$tag$rng,[$($args)*[$crate::rand_value!(@repeat$rng,[$ty]$len)]]$(,$($t)*)?)};(@$tag:ident$rng:expr,[$($args:tt)*][$($tt:tt)*]$(,$($t:tt)*)?)=>{$crate::rand_value!(@$tag$rng,[$($args)*[$crate::rand_value!(@repeat$rng,[$($tt)*])]]$(,$($t)*)?)};(@$tag:ident$rng:expr,[$($args:tt)*]$ty:expr)=>{$crate::rand_value!(@$tag$rng,[$($args)*[($rng).random($ty)]])};(@$tag:ident$rng:expr,[$($args:tt)*]$ty:expr,$($t:tt)*)=>{$crate::rand_value!(@$tag$rng,[$($args)*[($rng).random($ty)]],$($t)*)};(@$tag:ident$rng:expr,[$($args:tt)*],$($t:tt)*)=>{$crate::rand_value!(@$tag$rng,[$($args)*]$($t)*)};(@$tag:ident$rng:expr,[$($args:tt)*])=>{::std::compile_error!(::std::stringify!($($args)*))};($rng:expr,$($t:tt)*)=>{$crate::rand_value!(@inner$rng,[]$($t)*)}}#[macro_export]#[doc=" Declare random values using [`RandomSpec`]."]macro_rules!rand{(@assert$p:pat)=>{};(@assert$($p:tt)*)=>{::std::compile_error!(::std::concat!("expected pattern, found `",::std::stringify!($($p)*),"`"));};(@pat$rng:expr,[][])=>{};(@pat$rng:expr,[][],$($t:tt)*)=>{$crate::rand!(@pat$rng,[][]$($t)*)};(@pat$rng:expr,[$($p:tt)*][]$x:ident$($t:tt)*)=>{$crate::rand!(@pat$rng,[$($p)*$x][]$($t)*)};(@pat$rng:expr,[$($p:tt)*][]::$($t:tt)*)=>{$crate::rand!(@pat$rng,[$($p)*::][]$($t)*)};(@pat$rng:expr,[$($p:tt)*][]&$($t:tt)*)=>{$crate::rand!(@pat$rng,[$($p)*&][]$($t)*)};(@pat$rng:expr,[$($p:tt)*][]($($x:tt)*)$($t:tt)*)=>{$crate::rand!(@pat$rng,[$($p)*($($x)*)][]$($t)*)};(@pat$rng:expr,[$($p:tt)*][][$($x:tt)*]$($t:tt)*)=>{$crate::rand!(@pat$rng,[$($p)*[$($x)*]][]$($t)*)};(@pat$rng:expr,[$($p:tt)*][]{$($x:tt)*}$($t:tt)*)=>{$crate::rand!(@pat$rng,[$($p)*{$($x)*}][]$($t)*)};(@pat$rng:expr,[$($p:tt)*][]:$($t:tt)*)=>{$crate::rand!(@ty$rng,[$($p)*][]$($t)*)};(@ty$rng:expr,[$($p:tt)*][$($tt:tt)*]($($x:tt)*)$(,$($t:tt)*)?)=>{$crate::rand!(@let$rng,[$($p)*][$($tt)*($($x)*)]$(,$($t)*)?)};(@ty$rng:expr,[$($p:tt)*][$($tt:tt)*][$($x:tt)*]$(,$($t:tt)*)?)=>{$crate::rand!(@let$rng,[$($p)*][$($tt)*[$($x)*]]$(,$($t)*)?)};(@ty$rng:expr,[$($p:tt)*][$($tt:tt)*]$e:expr)=>{$crate::rand!(@let$rng,[$($p)*][$($tt)*$e])};(@ty$rng:expr,[$($p:tt)*][$($tt:tt)*]$e:expr,$($t:tt)*)=>{$crate::rand!(@let$rng,[$($p)*][$($tt)*$e],$($t)*)};(@ty$rng:expr,[$($p:tt)*][$($tt:tt)*]$e:tt)=>{$crate::rand!(@let$rng,[$($p)*][$($tt)*$e])};(@ty$rng:expr,[$($p:tt)*][$($tt:tt)*]$e:tt,$($t:tt)*)=>{$crate::rand!(@let$rng,[$($p)*][$($tt)*$e],$($t)*)};(@let$rng:expr,[$($p:tt)*][$($tt:tt)*]$($t:tt)*)=>{$crate::rand!{@assert$($p)*}let$($p)* =$crate::rand_value!($rng,$($tt)*);$crate::rand!(@pat$rng,[][]$($t)*)};($rng:expr)=>{};($rng:expr,$($t:tt)*)=>{$crate::rand!(@pat$rng,[][]$($t)*)};}} // 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;type Multiplicative:Monoid;#[doc=" additive identity: $0$"]fn zero()->Self::T{::unit()}fn is_zero(x:&Self::T)->bool where Self::T:PartialEq{*x==Self::zero()}#[doc=" multiplicative identity: $1$"]fn one()->Self::T{::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{::operate(x,y)}#[doc=" multiplicative operaion: $+$"]fn mul(x:&Self::T,y:&Self::T)->Self::T{::operate(x,y)}fn try_matrix_product(_a:&[Vec],_b:&[Vec])->Option>>{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){::operate_assign(x,y);}fn mul_assign(x:&mut Self::T,y:&Self::T){::operate_assign(x,y);}}pub trait Ring:SemiRing{#[doc=" additive inverse: $-$"]fn neg(x:&Self::T)->Self::T{::inverse(x)}#[doc=" additive right inversed operaion: $-$"]fn sub(x:&Self::T,y:&Self::T)->Self::T{::rinv_operate(x,y)}fn sub_assign(x:&mut Self::T,y:&Self::T){::rinv_operate_assign(x,y);}}implRing for R where R:SemiRing{}pub trait Field:Ring{#[doc=" multiplicative inverse: $-$"]fn inv(x:&Self::T)->Self::T{::inverse(x)}#[doc=" multiplicative right inversed operaion: $-$"]fn div(x:&Self::T,y:&Self::T)->Self::T{::rinv_operate(x,y)}fn div_assign(x:&mut Self::T,y:&Self::T){::rinv_operate_assign(x,y);}}implField for F where F:Ring{}pub trait DotProduct:Sized+Clone+Zero+Add+Mul{fn try_matrix_product(_a:&[Vec],_b:&[Vec])->Option>>{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);implDotProduct for Wrappingwhere Wrapping:Clone+Zero+Add>+Mul>{}#[doc=" $+,\\times$"]pub struct AddMulOperationwhere T:DotProduct+One{_marker:PhantomDataT>}implSemiRing for AddMulOperationwhere T:DotProduct+One{type T=T;type Additive=AdditiveOperation;type Multiplicative=MultiplicativeOperation;#[inline]fn try_matrix_product(a:&[Vec],b:&[Vec])->Option>>{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.}{Zero zero std::num::Wrapping,Self(T::zero())}{One one std::num::Wrapping,Self(T::one())});}