// codesnip-guard: FormalPowerSeries #[cfg_attr(any(),rustfmt::skip)]pub use self::formal_power_series::{FormalPowerSeries,FormalPowerSeriesCoefficient,FormalPowerSeriesCoefficientSqrt,Fps,Fps998244353};#[cfg_attr(any(),rustfmt::skip)]mod formal_power_series{use super::{Convolve998244353,ConvolveSteps,MInt,MIntConvert,MIntConvolve,MemorizedFactorial,NttReuse,One,PartialIgnoredOrd,Zero,gcd,montgomery::MInt998244353};use std::{fmt::{self,Debug},marker::PhantomData,ops::{Add,AddAssign,Div,DivAssign,Mul,MulAssign,Neg,Sub,SubAssign}};#[derive(Default)]pub struct FormalPowerSeries{pub data:Vec,_marker:PhantomData}implDebug for FormalPowerSerieswhere T:Debug{fn fmt(&self,f:&mut fmt::Formatter<'_>)->fmt::Result{Debug::fmt(&self.data,f)}}pub type Fps998244353=FormalPowerSeries;pub type Fps=FormalPowerSeries,MIntConvolve>;pub trait FormalPowerSeriesCoefficient:Sized+Clone+Zero+PartialEq+One+From+From+Add+Sub+Mul+Div+for<'r>Add<&'r Self,Output=Self>+for<'r>Sub<&'r Self,Output=Self>+for<'r>Mul<&'r Self,Output=Self>+for<'r>Div<&'r Self,Output=Self>+AddAssign+SubAssign+MulAssign+DivAssign+for<'r>AddAssign<&'r Self>+for<'r>SubAssign<&'r Self>+for<'r>MulAssign<&'r Self>+for<'r>DivAssign<&'r Self>+Neg{type Base:MIntConvert+MIntConvert;fn pow(self,exp:usize)->Self;fn signed_pow(self,exp:isize)->Self{if exp>=0{self.pow(exp as usize)}else{Self::one()/self.pow((-exp)as usize)}}fn memorized_factorial(n:usize)->MemorizedFactorial{MemorizedFactorial::new(n)}fn memorized_fact(mf:&MemorizedFactorial)->&[Self];fn memorized_inv_fact(mf:&MemorizedFactorial)->&[Self];fn memorized_inv(mf:&MemorizedFactorial,n:usize)->Self;}implFormalPowerSeriesCoefficient for MIntwhere M:MIntConvert+MIntConvert{type Base=M;fn pow(self,exp:usize)->Self{Self::pow(self,exp)}fn memorized_fact(mf:&MemorizedFactorial)->&[Self]{&mf.fact}fn memorized_inv_fact(mf:&MemorizedFactorial)->&[Self]{&mf.inv_fact}fn memorized_inv(mf:&MemorizedFactorial,n:usize)->Self{mf.inv(n)}}pub trait FormalPowerSeriesCoefficientSqrt:FormalPowerSeriesCoefficient{fn sqrt_coefficient(&self)->Option;}implFormalPowerSeriesCoefficientSqrt for MIntwhere M:MIntConvert+MIntConvert+MIntConvert{fn sqrt_coefficient(&self)->Option{self.sqrt()}}mod berlekamp_massey{use super::{FormalPowerSeries,FormalPowerSeriesCoefficient,NttReuse,One,Zero};use std::mem::{replace,swap};#[derive(Clone)]struct FpsMatrix{a00:FormalPowerSeries,a01:FormalPowerSeries,a10:FormalPowerSeries,a11:FormalPowerSeries}struct FrequencyMatrixwhere C:NttReuse{a00:C::F,a01:C::F,a10:C::F,a11:C::F}implClone for FrequencyMatrixwhere C:NttReuse,C::F:Clone{fn clone(&self)->Self{Self{a00:self.a00.clone(),a01:self.a01.clone(),a10:self.a10.clone(),a11:self.a11.clone()}}}implFpsMatrixwhere T:FormalPowerSeriesCoefficient,C:NttReuse>,C::F:Clone{fn identity()->Self{Self{a00:FormalPowerSeries::one(),a01:FormalPowerSeries::zero(),a10:FormalPowerSeries::zero(),a11:FormalPowerSeries::one()}}fn multiply_vector(&self,p:&FormalPowerSeries,q:&FormalPowerSeries)->(FormalPowerSeries,FormalPowerSeries){(add(&self.a00*p,&self.a01*q),add(&self.a10*p,&self.a11*q))}fn left_multiply_step(&mut self,quotient:&[T]){swap(&mut self.a00,&mut self.a10);swap(&mut self.a01,&mut self.a11);let quotient:FormalPowerSeries=FormalPowerSeries::from_vec(quotient.to_vec());self.a10=add(replace(&mut self.a10,FormalPowerSeries::zero()),quotient.clone()*&self.a00);self.a11=add(replace(&mut self.a11,FormalPowerSeries::zero()),quotient*&self.a01);}fn transform(&self,length:usize)->FrequencyMatrix{FrequencyMatrix{a00:reduced_transform(&self.a00,length),a01:reduced_transform(&self.a01,length),a10:reduced_transform(&self.a10,length),a11:reduced_transform(&self.a11,length)}}fn extend_transform(&self,frequency:FrequencyMatrix,length:usize)->FrequencyMatrix{fn extend(fps:&FormalPowerSeries,frequency:C::F,length:usize)->C::F where T:FormalPowerSeriesCoefficient,C:NttReuse>,C::F:Clone{if fps.length()<=length/2{C::ntt_doubling(frequency,false)}else{reduced_transform(fps,length)}}FrequencyMatrix{a00:extend(&self.a00,frequency.a00,length),a01:extend(&self.a01,frequency.a01,length),a10:extend(&self.a10,frequency.a10,length),a11:extend(&self.a11,frequency.a11,length)}}}implFrequencyMatrixwhere T:FormalPowerSeriesCoefficient,C:NttReuse>,C::F:Clone{fn product_sum(left_a:&C::F,right_a:&C::F,left_b:&C::F,right_b:&C::F)->C::F{let mut result=left_a.clone();C::multiply_prefix(&mut result,right_a);C::multiply_add(&mut result,left_b,right_b);result}fn multiply(&self,right:&Self)->Self{Self{a00:Self::product_sum(&self.a00,&right.a00,&self.a01,&right.a10),a01:Self::product_sum(&self.a00,&right.a01,&self.a01,&right.a11),a10:Self::product_sum(&self.a10,&right.a00,&self.a11,&right.a10),a11:Self::product_sum(&self.a10,&right.a01,&self.a11,&right.a11)}}fn apply(&self,p:&C::F,q:&C::F,length:usize)->(Vec,Vec){(C::inverse_transform_ntt(Self::product_sum(p,&self.a00,q,&self.a01),length),C::inverse_transform_ntt(Self::product_sum(p,&self.a10,q,&self.a11),length))}fn left_multiply_step(self,quotient:&FormalPowerSeries,length:usize)->Self{let negative_quotient=reduced_transform(&(-quotient),length);let mut a10=self.a00;C::multiply_add(&mut a10,&negative_quotient,&self.a10);let mut a11=self.a01;C::multiply_add(&mut a11,&negative_quotient,&self.a11);let result=Self{a00:self.a10,a01:self.a11,a10,a11};if C::MULTIPLE{result.inverse_transform(length).transform(length)}else{result}}fn inverse_transform(self,length:usize)->FpsMatrix{FpsMatrix{a00:FormalPowerSeries::from_vec(C::inverse_transform_ntt(self.a00,length)),a01:FormalPowerSeries::from_vec(C::inverse_transform_ntt(self.a01,length)),a10:FormalPowerSeries::from_vec(C::inverse_transform_ntt(self.a10,length)),a11:FormalPowerSeries::from_vec(C::inverse_transform_ntt(self.a11,length))}}}fn berlekamp_massey_naive(a:&[T],max_work:usize)->Option>where T:FormalPowerSeriesCoefficient{let n=a.len();let mut b=Vec::with_capacity(n+1);let mut c=Vec::with_capacity(n+1);let mut temporary=Vec::with_capacity(n+1);b.push(T::one());c.push(T::one());let mut y=T::one();let mut work=0usize;for k in 1..=n{let c_len=c.len();work=work.saturating_add(c_len);if work>max_work{return None;}let mut x=T::zero();for(c,a)in c.iter().zip(&a[k-c_len..]){x+=c.clone()*a.clone();}b.push(T::zero());let b_len=b.len();if x.is_zero(){continue;}let frequency=x.clone()/y.clone();if c_lenFormalPowerSerieswhere T:FormalPowerSeriesCoefficient,C:NttReuse>,C::F:Clone{pub fn berlekamp_massey(input:&[T])->Self{if input.last().is_none_or(|value|value.is_zero())&&input.iter().all(|value|value.is_zero()){return Self::one();}let max_work=if input.len()<=1536{usize::MAX}else{input.len().saturating_mul(2)};if let Some(recurrence)=berlekamp_massey_naive(input,max_work){return Self::from_vec(recurrence);}let n=input.len();let leading_zeros=input.iter().take_while(|value|value.is_zero()).count();let sequence=Self::from_vec(input.to_vec()).trimed();let mut modulus=Self::zeros(n+1);modulus[n]=T::one();let(matrix,_)=half_gcd(&modulus,&sequence,n/2,n.max(1).next_power_of_two());let(x,y)=matrix.multiply_vector(&modulus,&sequence);let mut recurrence=if y.length()==0{matrix.a01.clone()}else{matrix.a11.clone()};let recurrence_leading_zeros=recurrence.iter().take_while(|value|value.is_zero()).count();if recurrence_leading_zeros>0{let(division,_)=x.div_rem(y.clone());recurrence=add(recurrence*division,matrix.a01);}let inverse=T::one()/&recurrence[0];for value in recurrence.iter_mut(){*value*=&inverse;}let minimum_length=(leading_zeros+2).max(y.length()+1);if recurrence.length()(fps:&FormalPowerSeries)->isize{fps.length()as isize-1}fn add(left:FormalPowerSeries,right:FormalPowerSeries)->FormalPowerSerieswhere T:FormalPowerSeriesCoefficient{(left+right).trimed()}fn tail(fps:&FormalPowerSeries,start:isize)->FormalPowerSerieswhere T:FormalPowerSeriesCoefficient{let start=start.max(0)as usize;if start>=fps.length(){FormalPowerSeries::zero()}else{FormalPowerSeries::from_vec(fps.data[start..].to_vec())}}fn coefficient(fps:&FormalPowerSeries,index:isize)->T where T:FormalPowerSeriesCoefficient{if index<0{T::zero()}else{fps.coeff(index as usize)}}fn brute_force(mut p:FormalPowerSeries,mut q:FormalPowerSeries,k:usize)->FpsMatrixwhere T:FormalPowerSeriesCoefficient,C:NttReuse>,C::F:Clone{let threshold=degree(&p)-k as isize;let mut matrix=FpsMatrix::identity();while q.length()as isize>threshold{let q_degree=q.length()-1;let mut negative_quotient=vec![T::zero();p.length()-q.length()+1];let inverse=-T::one()/&q[q_degree];for i in(0..negative_quotient.len()).rev(){negative_quotient[i]=p[i+q_degree].clone()*&inverse;p[i+q_degree]=T::zero();for j in 0..q_degree{let value=negative_quotient[i].clone()*&q[j];p[i+j]+=&value;}}matrix.left_multiply_step(&negative_quotient);p.truncate(q_degree);p.trim_tail_zeros();swap(&mut p,&mut q);}matrix}fn reduced_transform(fps:&FormalPowerSeries,length:usize)->C::F where T:FormalPowerSeriesCoefficient,C:NttReuse>{let mut coefficients=vec![T::zero();length];for(i,value)in fps.iter().enumerate(){coefficients[i&(length-1)]+=value;}C::transform_ntt(coefficients,length)}fn transform_window(fps:&FormalPowerSeries,end:isize,length:usize)->C::F where T:FormalPowerSeriesCoefficient,C:NttReuse>{let start=end-length as isize;let coefficients=(start..end).map(|index|coefficient(fps,index)).collect();C::transform_ntt(coefficients,length)}fn half_gcd(p:&FormalPowerSeries,q:&FormalPowerSeries,k:usize,length:usize)->(FpsMatrix,FrequencyMatrix)where T:FormalPowerSeriesCoefficient,C:NttReuse>,C::F:Clone{let d=degree(p);if degree(q)=FormalPowerSeries::from_vec(q_reduced).trimed();let position=d-half as isize+degeneracy;let mut leading=T::zero();for i in 0..=position{leading+=coefficient(p,i)*coefficient(&matrix.a00,position-i)+coefficient(q,i)*coefficient(&matrix.a01,position-i);}p_reduced.push(leading);let mut p_reduced:FormalPowerSeries=FormalPowerSeries::from_vec(p_reduced);if degree(&q_reduced)<3*half as isize+degeneracy-k as isize{return(matrix,matrix_frequency);}let mut remaining=k as isize-degree(&matrix.a11);let mut top_product=matrix.a11.data.last().unwrap().clone();let mut product_degree=degree(&matrix.a11);if degeneracy>0{let skip=(2*half as isize+2*degeneracy-(d-half as isize+degeneracy)).max(0);let(division,remainder)=tail(&p_reduced,skip).div_rem(tail(&q_reduced,skip));remaining-=degree(&division);top_product*=-division.data.last().unwrap().clone();product_degree+=degree(&division);matrix_frequency=matrix_frequency.left_multiply_step(&division,length);swap(&mut p_reduced,&mut q_reduced);q_reduced=FormalPowerSeries::zeros(skip as usize);q_reduced.data.extend(remainder.data);}let start=3*half as isize+degeneracy-k as isize-remaining;let(right_matrix,right_frequency)=half_gcd(&tail(&p_reduced,start),&tail(&q_reduced,start),remaining as usize,length);let product_frequency=right_frequency.multiply(&matrix_frequency);let mut product=product_frequency.clone().inverse_transform(length);product.a00.truncate(k);product.a00.trim_tail_zeros();product.a01.truncate(k);product.a01.trim_tail_zeros();product.a10.truncate(k);product.a10.trim_tail_zeros();product_degree+=degree(&right_matrix.a11);if product_degree==length as isize{product.a11.resize(k+1);let highest=top_product*right_matrix.a11.data.last().unwrap();product.a11[k]=highest.clone();product.a11[0]-=highest;}product.a11.trim_tail_zeros();let product_frequency=if C::MULTIPLE{product.transform(length)}else{product_frequency};(product,product_frequency)}}mod formal_power_series_impls{use super::*;use std::{cmp::Reverse,collections::BinaryHeap,iter::repeat_with,iter::{FromIterator,once},marker::PhantomData,ops::{Index,IndexMut},slice::{Iter,IterMut}};implFormalPowerSeries{pub fn from_vec(data:Vec)->Self{Self{data,_marker:PhantomData}}pub fn length(&self)->usize{self.data.len()}pub fn truncate(&mut self,deg:usize){self.data.truncate(deg)}pub fn iter(&self)->Iter<'_,T>{self.data.iter()}pub fn iter_mut(&mut self)->IterMut<'_,T>{self.data.iter_mut()}}implClone for FormalPowerSerieswhere T:Clone{fn clone(&self)->Self{Self::from_vec(self.data.clone())}}implPartialEq for FormalPowerSerieswhere T:PartialEq{fn eq(&self,other:&Self)->bool{self.data.eq(&other.data)}}implEq for FormalPowerSerieswhere T:PartialEq{}implFormalPowerSerieswhere T:Zero{pub fn zeros(deg:usize)->Self{repeat_with(T::zero).take(deg).collect()}pub fn resize(&mut self,deg:usize){self.data.resize_with(deg,Zero::zero)}pub fn resized(mut self,deg:usize)->Self{self.resize(deg);self}pub fn reversed(mut self)->Self{self.data.reverse();self}}implFormalPowerSerieswhere T:Zero+Clone{pub fn coeff(&self,deg:usize)->T{self.data.get(deg).cloned().unwrap_or_else(T::zero)}}implFormalPowerSerieswhere T:Zero+PartialEq{pub fn trim_tail_zeros(&mut self){let mut len=self.length();while len>0{if self.data[len-1].is_zero(){len-=1;}else{break;}}self.truncate(len);}pub fn trimed(mut self)->Self{self.trim_tail_zeros();self}}implZero for FormalPowerSerieswhere T:PartialEq{fn zero()->Self{Self::from_vec(Vec::new())}}implOne for FormalPowerSerieswhere T:PartialEq+One{fn one()->Self{Self::from(T::one())}}implIntoIterator for FormalPowerSeries{type Item=T;type IntoIter=std::vec::IntoIter;fn into_iter(self)->Self::IntoIter{self.data.into_iter()}}impl<'a,T,C>IntoIterator for&'a FormalPowerSeries{type Item=&'a T;type IntoIter=Iter<'a,T>;fn into_iter(self)->Self::IntoIter{self.data.iter()}}impl<'a,T,C>IntoIterator for&'a mut FormalPowerSeries{type Item=&'a mut T;type IntoIter=IterMut<'a,T>;fn into_iter(self)->Self::IntoIter{self.data.iter_mut()}}implFromIteratorfor FormalPowerSeries{fn from_iter>(iter:I)->Self{Self::from_vec(iter.into_iter().collect())}}implIndexfor FormalPowerSeries{type Output=T;fn index(&self,index:usize)->&Self::Output{&self.data[index]}}implIndexMutfor FormalPowerSeries{fn index_mut(&mut self,index:usize)->&mut Self::Output{&mut self.data[index]}}implFromfor FormalPowerSeries{fn from(x:T)->Self{once(x).collect()}}implFrom>for FormalPowerSeries{fn from(data:Vec)->Self{Self::from_vec(data)}}implFormalPowerSerieswhere T:FormalPowerSeriesCoefficient{pub fn prefix_ref(&self,deg:usize)->Self{if degSelf{self.data.truncate(deg);self}pub fn even(mut self)->Self{let mut keep=false;self.data.retain(|_|{keep=!keep;keep});self}pub fn odd(mut self)->Self{let mut keep=true;self.data.retain(|_|{keep=!keep;keep});self}pub fn diff(mut self)->Self{let mut c=T::one();for i in 1..self.length(){self.data[i-1]=self.data[i].clone()*&c;c+=T::one();}self.data.pop();self}pub fn integral(mut self)->Self{let n=self.length();let mut fact=Vec::with_capacity(n+1);let mut c=T::one();fact.push(c.clone());for _ in 1..n{fact.push(fact.last().cloned().unwrap()*c.clone());c+=T::one();}let mut invf=T::one()/(fact.last().cloned().unwrap()*c.clone());self.data.push(T::zero());for i in(1..=n).rev(){self.data[i]=self.data[i-1].clone()*(invf.clone()*fact.pop().unwrap());invf*=c.clone();c-=T::one();}self.data[0]=T::zero();self}pub fn parity_inversion(mut self)->Self{self.iter_mut().skip(1).step_by(2).for_each(|x|*x=-x.clone());self}pub fn eval(&self,x:T)->T{self.iter().rev().fold(T::zero(),|sum,a|x.clone()*sum+a.clone())}}implFormalPowerSerieswhere T:FormalPowerSeriesCoefficient,C:ConvolveSteps>{#[inline]fn sparse_stride(&self,deg:usize,factor:usize)->Option{let limit=deg.next_power_of_two().trailing_zeros()as usize*factor;let mut count=0;let mut step=0;for(i,value)in self.iter().take(deg).enumerate(){if value.is_zero(){continue;}count+=1;if step!=1{step=gcd(step,i as u64);}if count>limit{return None;}}Some(step.max(1)as usize)}pub fn inv(&self,deg:usize)->Self{if deg==0{return Self::zero();}debug_assert!(!self[0].is_zero());if let Some(step)=self.sparse_stride(deg,6){let inv=T::one()/self[0].clone();let pos:Vec<_>=self.data.iter().take(deg).enumerate().skip(1).filter(|(_,x)|!x.is_zero()).map(|(i,x)|(i,-x.clone()*&inv)).collect();let mut f=Self::zeros(deg);f[0]=inv;for i in(pos.first().map_or(deg,|x|x.0)..deg).step_by(step){let mut tot=T::zero();for(j,coefficient)in&pos{if*j>i{break;}tot+=coefficient.clone()*&f[i-*j];}f[i]=tot;}return f;}let mut f=Self::from(T::one()/self[0].clone());f.data.reserve(deg.saturating_sub(1));let extend=|f:&mut Self,end|{for i in f.length()..end{let mut tot=T::zero();for j in 1..=i.min(self.length()-1){tot+=self[j].clone()*&f[i-j];}f.data.push(-tot*&f[0]);}};extend(&mut f,deg.min(32));let mut error=Vec::new();let mut i=f.length();while iSelf where C:NttReuse>,C::F:Clone{if deg==0{return Self::zero();}debug_assert!(self[0].is_zero());if let Some(step)=self.sparse_stride(deg,if deg<=256{16}else{8}){let diff=self.prefix_ref(deg).diff();let pos:Vec<_>=diff.data.iter().enumerate().filter_map(|(i,x)|if x.is_zero(){None}else{Some(i)}).collect();let mut f=Self::zeros(deg);f[0]=T::one();if pos.is_empty(){return f;}let mf=T::memorized_factorial(deg);for i in(pos.first().map_or(deg,|j|j+1)..deg).step_by(step){let mut tot=T::zero();for&j in&pos{if j>i-1{break;}tot+=f[i-1-j].clone()*&diff[j];}f[i]=tot*T::memorized_inv(&mf,i);}return f;}self.exp_or_pow(None,deg)}fn sum_products(f:&[C::F],g:&[C::F],len:usize)->Vecwhere C:NttReuse>,C::F:Clone{let chunk=C::max_product_sum_count(&f[0]);f.rchunks(chunk).zip(g.chunks(chunk)).map(|(f,g)|{let mut sum=f[f.len()-1].clone();C::multiply_prefix(&mut sum,&g[0]);for(f,g)in f.iter().rev().skip(1).zip(&g[1..]){C::multiply_add(&mut sum,f,g);}C::inverse_transform_ntt(sum,len)}).reduce(|mut sum,part|{for(sum,value)in sum.iter_mut().zip(part){*sum+=value;}sum}).unwrap()}fn exp_or_pow(&self,power:Option,deg:usize)->Self where C:NttReuse>,C::F:Clone{if deg==1{return Self::one();}let indices:Vec<_>=(0..=deg).map(T::from).collect();let modulus=>::mod_into();let mut inv=vec![T::zero();deg+1];inv[1]=T::one();for i in 2..=deg{inv[i]=-T::from(modulus/i)*&inv[modulus%i];}let block=deg.next_power_of_two()/16;let logarithm=if let Some(rhs)=&power{self.prefix_ref(block).log(block)*rhs}else{self.prefix_ref(block)};let(kernel,mut kernel_inverse,previous_inverse_fft)=logarithm.exp_newton(block,&indices,&inv);if power.is_some(){kernel_inverse=(self.prefix_ref(block)*&kernel).inv(block);}else{let mut error_fft=C::transform_ntt(kernel.data.clone(),block);C::multiply_prefix(&mut error_fft,&previous_inverse_fft);let error=C::inverse_transform_ntt(error_fft,block);let mut error_fft=C::transform_ntt(error.into_iter().skip(block/2).collect(),block);C::multiply_prefix(&mut error_fft,&previous_inverse_fft);let error=C::inverse_transform_ntt(error_fft,block/2);kernel_inverse.data.extend(error.into_iter().take(block/2).map(Neg::neg));}let kernel_data=kernel.data;let kernel_inverse_data=kernel_inverse.data;let kernel_inverse=C::transform(kernel_inverse_data,block*2);let kernel=C::transform(kernel_data.clone(),block*2);let blocks=deg.div_ceil(block);let mut derivative_ffts=Vec::with_capacity(blocks-1);let mut polynomial_ffts=Vec::with_capacity(if power.is_some(){blocks-1}else{0});for q in 1..blocks{let mut values=Self::zeros(block*2);for(i,values)in values.data.chunks_mut(block).enumerate(){let start=(q-i)*block;for(value,x)in values.iter_mut().zip(self.iter().skip(start).take(deg-start)){*value=x.clone();}}if power.is_some(){polynomial_ffts.push(C::transform_ntt(values.data.clone(),block*2));}for(i,values)in values.data.chunks_mut(block).enumerate(){let start=(q-i)*block;for(value,index)in values.iter_mut().zip(&indices[start..]){*value*=index;}}derivative_ffts.push(C::transform_ntt(values.data,block*2));}let mut result=kernel_data.clone();result.reserve(deg-block);let mut result_ffts=Vec::with_capacity(blocks-1);for q in 1..blocks{result_ffts.push(C::transform_ntt(result[(q-1)*block..q*block].to_vec(),block*2));let mut values=Self::sum_products(&derivative_ffts[..q],&result_ffts,block);if let Some(rhs)=&power{let product=Self::sum_products(&polynomial_ffts[..q],&result_ffts,block);let factor=rhs.clone()+T::one();for(i,value)in values.iter_mut().take(deg-q*block).enumerate(){*value=value.clone()*&factor-product[i].clone()*&indices[q*block+i];}}let mut values=C::transform(values,block*2);C::multiply(&mut values,&kernel_inverse);let mut values=C::inverse_transform(values,block*2);values.truncate(block);let len=block.min(deg-q*block);for(i,value)in values.iter_mut().take(len).enumerate(){*value*=&inv[q*block+i];}values[len..].fill(T::zero());let mut values=C::transform(values,block*2);C::multiply(&mut values,&kernel);let mut values=C::inverse_transform(values,block*2);values.truncate(len);result.extend(values);}Self::from_vec(result)}fn exp_newton(&self,deg:usize,indices:&[T],inv:&[T])->(Self,Self,C::F)where C:NttReuse>,C::F:Clone{if deg==1{let one=Self::one();return(one.clone(),one.clone(),C::transform_ntt(one.data,1));}let mut f=Self::from_vec(vec![T::one(),self.coeff(1)]);let mut inverse=Self::one();let mut inverse_fft=C::transform_ntt(inverse.data.clone(),2);let mut m=2;while mSelf{if deg==0{return Self::zero();}debug_assert!(!self[0].is_zero());if deg==1{return Self::zeros(1);}if let Some(step)=self.sparse_stride(deg,2){let pos:Vec<_>=self.iter().take(deg).enumerate().skip(1).filter_map(|(i,x)|(!x.is_zero()).then_some(i)).collect();let mut derivative=Self::zeros(deg);let inverse=T::one()/self[0].clone();for i in(pos.first().copied().unwrap_or(deg)..deg).step_by(step){let mut value=self.coeff(i)*T::from(i);for&j in&pos{if j>=i{break;}value-=self[j].clone()*&derivative[i-j];}derivative[i]=value*&inverse;}if pos.is_empty(){return derivative;}derivative.data.remove(0);return derivative.integral();}let n=deg-1;if n<=64{return(self.inv(deg)*self.prefix_ref(deg).diff()).prefix(n).integral();}let half=n.next_power_of_two()/2;let derivative=self.prefix_ref(deg).diff();let inverse=C::transform(self.inv(half).data,half*2);let mut quotient=C::transform(derivative.prefix_ref(half).data,half*2);C::multiply(&mut quotient,&inverse);let mut result=C::inverse_transform(quotient,half*2);result.truncate(half);if n-half<=4{let inverse=T::one()/self[0].clone();for i in half..n{let mut value=derivative.coeff(i);for j in 1..=i.min(self.length()-1){value-=self[j].clone()*&result[i-j];}result.push(value*&inverse);}return Self::from_vec(result).integral();}let quotient=C::transform(result.clone(),half*2);let mut error=C::transform(self.prefix_ref(n).data,half*2);C::multiply(&mut error,"ient);let mut error=C::inverse_transform(error,half*2);for i in 0..n-half{error[i]=derivative.coeff(half+i)-&error[half+i];}error.truncate(n-half);let mut error=C::transform(error,half*2);C::multiply(&mut error,&inverse);let error=C::inverse_transform(error,half*2);result.extend(error.into_iter().take(n-half));Self::from_vec(result).integral()}pub fn pow(&self,rhs:usize,deg:usize)->Self where C:NttReuse>,C::F:Clone{if rhs==0{return Self::from_vec(once(T::one()).chain(repeat_with(T::zero)).take(deg).collect());}if rhs==1{return self.prefix_ref(deg).resized(deg);}if let Some(k)=self.iter().take(deg.div_ceil(rhs)).position(|x|!x.is_zero()){let deg=deg-k*rhs;let x0=self[k].clone();let mut f=(self.prefix_ref(k+deg)>>k)/&x0;if let Some(step)=f.sparse_stride(deg,12){f=f.pow_sparse1(T::from(rhs),deg,step);}else if rhs<=4{let squared=(&f*&f).prefix(deg);f=match rhs{2=>squared,3=>(squared*f).prefix(deg),_=>(&squared*&squared).prefix(deg),}.resized(deg);}else{f=f.exp_or_pow(Some(T::from(rhs)),deg);}f*=x0.pow(rhs);f<<=k*rhs;f}else{Self::zeros(deg)}}fn pow_sparse1(&self,rhs:T,deg:usize,step:usize)->Self{debug_assert!(!self[0].is_zero());let mut pos:Vec<_>=self.data.iter().take(deg).enumerate().skip(1).filter(|(_,x)|!x.is_zero()).map(|(i,x)|(i,x.clone(),T::from(i)*&rhs*x)).collect();let mut f=Self::zeros(deg);f[0]=T::one();if pos.is_empty(){return f;}let mf=T::memorized_factorial(deg);for(_,coefficient,_)in&mut pos{*coefficient*=T::from(step);}for i in(pos.first().map_or(deg,|x|x.0)..deg).step_by(step){let mut tot=T::zero();for(j,coefficient,weight)in&mut pos{if*j>i{break;}tot+=weight.clone()*&f[i-*j];*weight-=&*coefficient;}f[i]=tot*T::memorized_inv(&mf,i);}f}fn sparse_fold(&self,sparse:impl IntoIterator,deg:usize)->T{sparse.into_iter().take_while(|&(i,_)|i<=deg).fold(T::zero(),|sum,(i,x)|sum+x*self.coeff(deg-i))}#[doc=" solve: $X(QF)'=\\alpha P'(QF)+\\beta P(Q'F)$ in $O(deg * max(nz(P), nz(Q), nz(X)))$"]pub fn solve_sparse_differential2(p:&Self,q:&Self,x:&Self,alpha:T,beta:T,deg:usize)->Self{if deg==0{return Self::zero();}let collect_sparse=|p:&Self|->Vec<(usize,T)>{p.iter().enumerate().filter(|&(_,x)|!x.is_zero()).map(|(i,x)|(i,x.clone())).collect()};assert!(q.coeff(0).is_one());assert!(x.coeff(0).is_one());let p=collect_sparse(p);let q=collect_sparse(q);let x=collect_sparse(x);let diff=|p:&[(usize,T)]|->Vec<(usize,T)>{p.iter().filter(|&&(i,_)|i>0).map(|&(i,ref x)|(i-1,x.clone()*T::from(i))).collect()};let dp=diff(&p);let dq=diff(&q);let mf=T::memorized_factorial(deg);let mut f=Self::zeros(deg);let mut qf=Self::zeros(deg);let mut dq_f=Self::zeros(deg);let mut d_qf=Self::zeros(deg);f[0]=T::one();for i in 0..deg-1{qf[i]=f.sparse_fold(q.iter().cloned(),i);dq_f[i]=f.sparse_fold(dq.iter().cloned(),i);let dp_qf_i=qf.sparse_fold(dp.iter().cloned(),i);let p_dq_f_i=dq_f.sparse_fold(p.iter().cloned(),i);let x_d_qf_i=d_qf.sparse_fold(x.iter().map(|&(i,ref x)|(i,x.clone()-T::from((i==0)as usize))),i);d_qf[i]=alpha.clone()*dp_qf_i+beta.clone()*p_dq_f_i-x_d_qf_i;let mut f_ip1=d_qf[i].clone();for&(j,ref q)in q.iter().take_while(|&&(j,_)|j<=i){if j>0{f_ip1-=q.clone()*&f[i-(j-1)]*T::from(i-(j-1));}}f[i+1]=f_ip1*T::memorized_inv(&mf,i+1);}f}#[doc=" P^exp_p * Q^exp_q"]pub fn mul_of_pow_sparse(&self,q:&Self,exp_p:isize,exp_q:isize,deg:usize)->Self{if deg==0{return Self::zero();}if exp_p==0&&exp_q==0{return Self::from_vec(once(T::one()).chain(repeat_with(T::zero)).take(deg).collect());}if exp_p!=0&&self.iter().all(|x|x.is_zero()){assert!(exp_p>0);return Self::zeros(deg);}if exp_q!=0&&q.iter().all(|x|x.is_zero()){assert!(exp_q>0);return Self::zeros(deg);}let normalize=|f:&Self,exp:isize|{if exp==0{return(0usize,T::one(),Self::from_vec(vec![T::one()]));}let k=f.iter().position(|value|!value.is_zero()).unwrap();assert!(exp>=0||k==0,"Negative exponent with zero constant term");let c=f[k].clone();let f=(f.clone()>>k)/&c;(k,c,f)};let(sp,cp,mut p)=normalize(self,exp_p);let(sq,cq,mut q)=normalize(q,exp_q);let shift=exp_p.saturating_mul(sp as _).saturating_add(exp_q.saturating_mul(sq as _))as usize;if shift>=deg{return Self::zeros(deg);}p.truncate(deg-shift);q.truncate(deg-shift);let mut f=Self::solve_sparse_differential2(&p,&q,&p,T::from(exp_p),T::from(exp_q),deg-shift);f*=cp.signed_pow(exp_p)*cq.signed_pow(exp_q);if shift>0{f<<=shift;}f.prefix(deg)}#[doc=" exp(P/Q)"]pub fn exp_of_div_sparse(&self,q:&Self,deg:usize)->Self{if deg==0{return Self::zero();}let shift_q=q.iter().position(|value|!value.is_zero()).expect("Zero denominator");let shift_p=self.iter().position(|value|!value.is_zero()).unwrap_or(!0);assert!(shift_p>shift_q);let mut p=self>>shift_q;let mut q=q>>shift_q;assert!(!q.coeff(0).is_zero());let c=q[0].clone();p/=c.clone();q/=c;Self::solve_sparse_differential2(&p,&q,&q,T::one(),-T::one(),deg)}}implFormalPowerSerieswhere T:FormalPowerSeriesCoefficientSqrt,C:ConvolveSteps>{pub fn sqrt(&self,deg:usize)->Option{if self[0].is_zero(){if let Some(k)=self.iter().position(|x|!x.is_zero()){if k%2!=0{return None;}else if deg>k/2{return Some((self>>k).sqrt(deg-k/2)?<<(k/2));}}}else{let s=self[0].sqrt_coefficient()?;if deg<=1{return Some(Self::from(s).prefix(deg));}if let Some(step)=self.sparse_stride(deg,4){let t=self[0].clone();let mut f=self.prefix_ref(deg)/t;f=f.pow_sparse1(T::one()/T::from(2usize),deg,step);f*=s;return Some(f);}let mut f=Self::from(s);let inv2=T::one()/(T::one()+T::one());let inv2s=inv2.clone()/&f[0];let extend=|f:&mut Self,end|{for i in f.length()..end{let mut value=self.coeff(i);for j in 1..i{value-=f[j].clone()*&f[i-j];}f.data.push(value*&inv2s);}};extend(&mut f,deg.min(32));f.truncate(deg);if f.length()==deg{return Some(f);}let mut inverse=f.inv(f.length());let mut i=f.length();while i>i}else{let square=C::square(f.data.clone(),i);Self::from_vec(square.into_iter().take(len-i).enumerate().map(|(j,value)|self.coeff(i+j)+self.coeff(j)-value).collect())};let mut error_fft=C::transform(error.data,i*2);C::multiply(&mut error_fft,&factor);let delta=C::inverse_transform(error_fft,i*2);f.data.extend(delta.into_iter().take(len-i).map(|x|x*&inv2));if i*2+4FormalPowerSerieswhere T:FormalPowerSeriesCoefficient,C:ConvolveSteps>{pub fn count_subset_sum(&self,deg:usize,mut inverse:F)->Self where F:FnMut(usize)->T,C:NttReuse>,C::F:Clone{let n=self.length();let mut f=Self::zeros(n);for i in 1..n{if!self[i].is_zero(){for(j,d)in(0..n).step_by(i).enumerate().skip(1){if j&1!=0{f[d]+=self[i].clone()*&inverse(j);}else{f[d]-=self[i].clone()*&inverse(j);}}}}f.exp(deg)}pub fn count_multiset_sum(&self,deg:usize,mut inverse:F)->Self where F:FnMut(usize)->T,C:NttReuse>,C::F:Clone{let n=self.length();let mut f=Self::zeros(n);for i in 1..n{if!self[i].is_zero(){for(j,d)in(0..n).step_by(i).enumerate().skip(1){f[d]+=self[i].clone()*&inverse(j);}}}f.exp(deg)}#[doc=" [x^n] P(x) / Q(x)"]pub fn bostan_mori(mut self,mut rhs:Self,mut n:usize)->T where C:NttReuse>{let mut res=T::zero();rhs.trim_tail_zeros();if self.length()>=rhs.length(){let r=&self/&rhs;if n0{let t=C::even_mul_normal_neg(&q,&q);p=if n.is_multiple_of(2){C::even_mul_normal_neg(&p,&q)}else{C::odd_mul_normal_neg(&p,&q)};q=t;n/=2;if n!=0{if nSelf{let d=self.length()-1;if n==0{return(Self::one()<<(d-1))/self[0].clone();}let q=self;let mq=q.clone().parity_inversion();let w=(q*&mq).even().bostan_mori_msb(n/2);let mut s=Self::zeros(w.length()*2-(n%2));for(i,x)in w.iter().enumerate(){s[i*2+(1-n%2)]=x.clone();}let len=2*d+1;let ts=C::transform(s.prefix(len).data,len);mq.reversed().middle_product(&ts,len).prefix(d+1)}#[doc=" x^n mod self"]pub fn pow_mod(self,n:usize)->Self{let d=self.length()-1;let q=self.reversed();let u=q.clone().bostan_mori_msb(n);let mut f=(u*q).prefix(d).reversed();f.trim_tail_zeros();f}fn middle_product(self,other:&C::F,deg:usize)->Self{let n=self.length();let mut s=C::transform(self.reversed().data,deg);C::multiply(&mut s,other);Self::from_vec((C::inverse_transform(s,deg))[n-1..].to_vec())}pub fn multipoint_evaluation(self,points:&[T])->Vecwhere C:NttReuse>,C::F:Clone{let n=points.len();if n<=32||self.length()<=32{return points.iter().map(|p|self.eval(p.clone())).collect();}let size=n.next_power_of_two();let block=16;let leaves=size/block;let mut subproduct_tree=Vec::with_capacity(leaves*2);subproduct_tree.resize_with(leaves*2,| |None);let mut leaf_products=Vec::with_capacity(leaves);for i in 0..leaves{let mut product=vec![T::one()];for j in 0..block{let x=points.get(i*block+j).cloned().unwrap_or_else(T::zero);product.push(T::one());for k in(1..=j).rev(){product[k]=product[k-1].clone()-x.clone()*&product[k];}product[0]*=-x;}subproduct_tree[leaves+i]=Some(C::transform_ntt(product.clone(),block*2));leaf_products.push(product);}for i in(1..leaves).rev(){let mut product=subproduct_tree[i*2].as_ref().unwrap().clone();C::multiply_prefix(&mut product,subproduct_tree[i*2+1].as_ref().unwrap());if i>1{product=C::ntt_doubling(product,true);}subproduct_tree[i]=Some(product);}let mut product=C::inverse_transform_ntt(subproduct_tree[1].take().unwrap(),size);product[0]-=T::one();product.push(T::one());let mut uptree_t=Vec::with_capacity(leaves*2);uptree_t.resize_with(1,Zero::zero);let m=self.length();let v=Self::from_vec(product).reversed().resized(m);let s=C::transform(self.data,m*2);uptree_t.push(v.inv(m).middle_product(&s,m*2).resized(size));for i in 1..leaves{let degree=uptree_t[i].length();let spectrum=C::transform_ntt(std::mem::take(&mut uptree_t[i].data),degree);let left=subproduct_tree[i*2].take().unwrap();let right=subproduct_tree[i*2+1].take().unwrap();let mut child=spectrum.clone();C::multiply_prefix(&mut child,&right);let mut child=C::inverse_transform_ntt(child,degree);child.drain(..degree/2);uptree_t.push(Self::from_vec(child));let mut child=spectrum;C::multiply_prefix(&mut child,&left);let mut child=C::inverse_transform_ntt(child,degree);child.drain(..degree/2);uptree_t.push(Self::from_vec(child));}let mut result=Vec::with_capacity(n);for((values,product),points)in uptree_t[leaves..].iter().zip(leaf_products).zip(points.chunks(block)){let mut remainder=Self::zeros(block);for(j,value)in values.iter().enumerate(){for(r,p)in remainder.data[..=j].iter_mut().zip(&product[block-j..]){*r+=value.clone()*p;}}result.extend(points.iter().map(|p|remainder.eval(p.clone())));}result}pub fn product_all(iter:I,deg:usize)->Self where I:IntoIterator{let mut heap:BinaryHeap<_>=iter.into_iter().map(|f|PartialIgnoredOrd(Reverse(f.length()),f)).collect();while let Some(PartialIgnoredOrd(_,x))=heap.pop(){if let Some(PartialIgnoredOrd(_,y))=heap.pop(){let z=(x*y).prefix(deg);heap.push(PartialIgnoredOrd(Reverse(z.length()),z));}else{return x;}}Self::one()}pub fn sum_all_rational(iter:I,deg:usize)->(Self,Self)where I:IntoIterator{let mut heap:BinaryHeap<_>=iter.into_iter().map(|(f,g)|PartialIgnoredOrd(Reverse(f.length().max(g.length())),(f,g))).collect();while let Some(PartialIgnoredOrd(_,(xa,xb)))=heap.pop(){if let Some(PartialIgnoredOrd(_,(ya,yb)))=heap.pop(){let zb=(&xb*&yb).prefix(deg);let za=(xa*yb+ya*xb).prefix(deg);heap.push(PartialIgnoredOrd(Reverse(za.length().max(zb.length())),(za,zb)));}else{return(xa,xb);}}(Self::zero(),Self::one())}pub fn kth_term_of_linearly_recurrence(self,a:Vec,k:usize)->T where C:NttReuse>{if let Some(x)=a.get(k){return x.clone();}let p=(Self::from_vec(a).prefix(self.length()-1)*&self).prefix(self.length()-1);p.bostan_mori(self,k)}pub fn kth_term(a:Vec,k:usize)->T where C:NttReuse>,C::F:Clone{if let Some(x)=a.get(k){return x.clone();}Self::berlekamp_massey(&a).kth_term_of_linearly_recurrence(a,k)}#[doc=" sum_i a_i exp(b_i x)"]pub fn linear_sum_of_exp(iter:I,deg:usize,mut inv_fact:F)->Self where I:IntoIterator,F:FnMut(usize)->T{let(p,q)=Self::sum_all_rational(iter.into_iter().map(|(a,b)|(Self::from_vec(vec![a]),Self::from_vec(vec![T::one(),-b]))),deg);let mut f=(p*q.inv(deg)).prefix(deg);for i in 0..f.length(){f[i]*=inv_fact(i);}f}#[doc=" sum_i (a_i x)^j"]pub fn sum_of_powers(iter:I,deg:usize)->Self where I:IntoIterator{let mut n=T::zero();let prod=Self::product_all(iter.into_iter().map(|a|{n+=T::one();Self::from_vec(vec![T::one(),-a])}),deg);(-prod.log(deg).diff()<<1)+Self::from_vec(vec![n])}pub fn power_projection(&self,w:&[T],m:usize)->Self where C:NttReuse>{if w.is_empty(){return Self::zeros(m);}if m<=1{return Self::from_vec(vec![w[0].clone();m]);}let n0=w.len();let mut n=n0.next_power_of_two();let mut f=self.prefix_ref(n);f.resize(n);let base=n*2;let mut p_flat=vec![T::zero();base];for(i,wi)in w.iter().enumerate(){p_flat[n-1-i]=wi.clone();}let mut q_flat=vec![T::zero();base*2];q_flat[0]=T::one();let q_offset=base;for(i,fi)in f.iter().enumerate(){q_flat[q_offset+i]=-fi.clone();}let mut py=1usize;let mut qy=2usize;let y_limit=m;while n>1{let(mut p,mut q)=C::power_projection_step(p_flat,q_flat,n,py,qy);let new_py=(py+qy-1).min(y_limit);let new_qy=(qy+qy-1).min(y_limit);p.resize_with(n*new_py,T::zero);q.resize_with(n*new_qy,T::zero);let n2=n/2;for row in p.chunks_exact_mut(n){row[n2..].fill_with(T::zero);}for row in q.chunks_exact_mut(n){row[n2..].fill_with(T::zero);}p_flat=p;q_flat=q;py=new_py;qy=new_qy;n=n2;}let base=2;let mut p_y=Vec::with_capacity(py);for y in 0..py{p_y.push(p_flat[base*y].clone());}let mut q_y=Vec::with_capacity(qy);for y in 0..qy{q_y.push(q_flat[base*y].clone());}(Self::from_vec(p_y)*Self::from_vec(q_y).inv(m)).prefix(m)}pub fn compositional_inverse(&self,deg:usize)->Self where C:NttReuse>,C::F:Clone{if deg==0{return Self::zero();}if deg==1{return Self::from_vec(vec![T::zero()]);}debug_assert!(self[0].is_zero());debug_assert!(!self[1].is_zero());let mut f=self.prefix_ref(deg);f.resize(deg);let c=f[1].clone();f/=c.clone();let mut w=vec![T::zero();deg];w[deg-1]=T::one();let s=f.power_projection(&w,deg);let n=deg-1;let n_t=T::from(n);let mut h=vec![T::zero();n];for i in 1..=n{h[n-i]=s[i].clone()*&n_t/T::from(i);}let h_fps=Self::from_vec(h);let inv_n=T::one()/n_t;let mut t=h_fps.log(n);t*=-inv_n;let g_over_x=t.exp(n);let mut g=(g_over_x<<1).prefix(deg);let inv_c=T::one()/c;let mut pow=T::one();for coef in g.iter_mut(){*coef*=pow.clone();pow*=inv_c.clone();}g}#[doc=" f(x) <- f(x + a)"]pub fn taylor_shift(mut self,a:T)->Self{let f=T::memorized_factorial(self.length());let n=self.length();for(i,coef)in self.data.iter_mut().enumerate(){*coef*=T::memorized_fact(&f)[i].clone();}self.data.reverse();let mut b=a.clone();let mut g=Self::from_vec(T::memorized_inv_fact(&f)[..n].to_vec());for i in 1..n{g[i]*=b.clone();b*=a.clone();}self*=g;self.truncate(n);self.data.reverse();for(i,coef)in self.data.iter_mut().enumerate(){*coef*=T::memorized_inv_fact(&f)[i].clone();}self}}}mod formal_power_series_nums{use super::*;use std::{mem::take,ops::{Add,AddAssign,Div,DivAssign,Mul,MulAssign,Neg,Rem,RemAssign,Shl,ShlAssign,Shr,ShrAssign,Sub,SubAssign}};implAddAssignfor FormalPowerSerieswhere T:FormalPowerSeriesCoefficient{fn add_assign(&mut self,rhs:T){if self.length()==0{self.data.push(T::zero());}self.data[0].add_assign(rhs);}}implSubAssignfor FormalPowerSerieswhere T:FormalPowerSeriesCoefficient{fn sub_assign(&mut self,rhs:T){if self.length()==0{self.data.push(T::zero());}self.data[0].sub_assign(rhs);self.trim_tail_zeros();}}implMulAssignfor FormalPowerSerieswhere T:FormalPowerSeriesCoefficient{fn mul_assign(&mut self,rhs:T){for x in self.iter_mut(){x.mul_assign(&rhs);}}}implDivAssignfor FormalPowerSerieswhere T:FormalPowerSeriesCoefficient{fn div_assign(&mut self,rhs:T){let rinv=T::one()/rhs;for x in self.iter_mut(){x.mul_assign(&rinv);}}}macro_rules!impl_fps_single_binop{($imp:ident,$method:ident,$imp_assign:ident,$method_assign:ident)=>{impl$imp_assign<&T>for FormalPowerSerieswhere T:FormalPowerSeriesCoefficient,{fn$method_assign(&mut self,rhs:&T){$imp_assign::$method_assign(self,rhs.clone());}}impl$impfor FormalPowerSerieswhere T:FormalPowerSeriesCoefficient,{type Output=Self;fn$method(mut self,rhs:T)->Self::Output{$imp_assign::$method_assign(&mut self,rhs);self}}impl$imp<&T>for FormalPowerSerieswhere T:FormalPowerSeriesCoefficient,{type Output=Self;fn$method(mut self,rhs:&T)->Self::Output{$imp_assign::$method_assign(&mut self,rhs);self}}impl$impfor&FormalPowerSerieswhere T:FormalPowerSeriesCoefficient,{type Output=FormalPowerSeries;fn$method(self,rhs:T)->Self::Output{$imp::$method(self.clone(),rhs)}}impl$imp<&T>for&FormalPowerSerieswhere T:FormalPowerSeriesCoefficient,{type Output=FormalPowerSeries;fn$method(self,rhs:&T)->Self::Output{$imp::$method(self.clone(),rhs)}}};}impl_fps_single_binop!(Add,add,AddAssign,add_assign);impl_fps_single_binop!(Sub,sub,SubAssign,sub_assign);impl_fps_single_binop!(Mul,mul,MulAssign,mul_assign);impl_fps_single_binop!(Div,div,DivAssign,div_assign);implAddAssign<&Self>for FormalPowerSerieswhere T:FormalPowerSeriesCoefficient{fn add_assign(&mut self,rhs:&Self){if self.length()SubAssign<&Self>for FormalPowerSerieswhere T:FormalPowerSeriesCoefficient{fn sub_assign(&mut self,rhs:&Self){if self.length(){impl$imp_assign for FormalPowerSerieswhere T:FormalPowerSeriesCoefficient,{fn$method_assign(&mut self,rhs:Self){$imp_assign::$method_assign(self,&rhs);}}impl$imp for FormalPowerSerieswhere T:FormalPowerSeriesCoefficient,{type Output=Self;fn$method(mut self,rhs:Self)->Self::Output{$imp_assign::$method_assign(&mut self,&rhs);self}}impl$imp<&FormalPowerSeries>for FormalPowerSerieswhere T:FormalPowerSeriesCoefficient,{type Output=Self;fn$method(mut self,rhs:&FormalPowerSeries)->Self::Output{$imp_assign::$method_assign(&mut self,rhs);self}}impl$imp>for&FormalPowerSerieswhere T:FormalPowerSeriesCoefficient,{type Output=FormalPowerSeries;fn$method(self,rhs:FormalPowerSeries)->Self::Output{let mut self_=self.clone();$imp_assign::$method_assign(&mut self_,&rhs);self_}}impl$imp<&FormalPowerSeries>for&FormalPowerSerieswhere T:FormalPowerSeriesCoefficient,{type Output=FormalPowerSeries;fn$method(self,rhs:&FormalPowerSeries)->Self::Output{let mut self_=self.clone();$imp_assign::$method_assign(&mut self_,rhs);self_}}};}impl_fps_binop_addsub!(Add,add,AddAssign,add_assign);impl_fps_binop_addsub!(Sub,sub,SubAssign,sub_assign);implMul for FormalPowerSerieswhere C:ConvolveSteps>{type Output=Self;fn mul(self,rhs:Self)->Self::Output{Self::from_vec(C::convolve(self.data,rhs.data))}}implDiv for FormalPowerSerieswhere T:FormalPowerSeriesCoefficient,C:ConvolveSteps>{type Output=Self;fn div(mut self,mut rhs:Self)->Self::Output{self.trim_tail_zeros();rhs.trim_tail_zeros();if self.length()Rem for FormalPowerSerieswhere T:FormalPowerSeriesCoefficient,C:ConvolveSteps>{type Output=Self;fn rem(self,rhs:Self)->Self::Output{let mut rem=self.clone()-self/rhs.clone()*rhs;rem.trim_tail_zeros();rem}}implFormalPowerSerieswhere T:FormalPowerSeriesCoefficient,C:ConvolveSteps>{pub fn div_rem(self,rhs:Self)->(Self,Self){let div=self.clone()/rhs.clone();let mut rem=self-div.clone()*rhs;rem.trim_tail_zeros();(div,rem)}}macro_rules!impl_fps_binop_conv{($imp:ident,$method:ident,$imp_assign:ident,$method_assign:ident)=>{impl$imp_assign for FormalPowerSerieswhere T:FormalPowerSeriesCoefficient,C:ConvolveSteps>,{fn$method_assign(&mut self,rhs:Self){*self=$imp::$method(Self::from_vec(take(&mut self.data)),rhs);}}impl$imp_assign<&Self>for FormalPowerSerieswhere T:FormalPowerSeriesCoefficient,C:ConvolveSteps>,{fn$method_assign(&mut self,rhs:&Self){$imp_assign::$method_assign(self,rhs.clone());}}impl$imp<&FormalPowerSeries>for FormalPowerSerieswhere T:FormalPowerSeriesCoefficient,C:ConvolveSteps>,{type Output=Self;fn$method(self,rhs:&FormalPowerSeries)->Self::Output{$imp::$method(self,rhs.clone())}}impl$imp>for&FormalPowerSerieswhere T:FormalPowerSeriesCoefficient,C:ConvolveSteps>,{type Output=FormalPowerSeries;fn$method(self,rhs:FormalPowerSeries)->Self::Output{$imp::$method(self.clone(),rhs)}}impl$imp<&FormalPowerSeries>for&FormalPowerSerieswhere T:FormalPowerSeriesCoefficient,C:ConvolveSteps>,{type Output=FormalPowerSeries;fn$method(self,rhs:&FormalPowerSeries)->Self::Output{$imp::$method(self.clone(),rhs.clone())}}};}impl_fps_binop_conv!(Mul,mul,MulAssign,mul_assign);impl_fps_binop_conv!(Div,div,DivAssign,div_assign);impl_fps_binop_conv!(Rem,rem,RemAssign,rem_assign);implNeg for FormalPowerSerieswhere T:FormalPowerSeriesCoefficient{type Output=Self;fn neg(mut self)->Self::Output{for x in self.iter_mut(){*x=-x.clone();}self}}implNeg for&FormalPowerSerieswhere T:FormalPowerSeriesCoefficient{type Output=FormalPowerSeries;fn neg(self)->Self::Output{self.clone().neg()}}implShrAssignfor FormalPowerSerieswhere T:FormalPowerSeriesCoefficient{fn shr_assign(&mut self,rhs:usize){if self.length()<=rhs{*self=Self::zero();}else{for i in rhs..self.length(){self[i-rhs]=self[i].clone();}self.truncate(self.length()-rhs);}}}implShlAssignfor FormalPowerSerieswhere T:FormalPowerSeriesCoefficient{fn shl_assign(&mut self,rhs:usize){let n=self.length();self.resize(n+rhs);for i in(0..n).rev(){self[i+rhs]=self[i].clone();}for i in 0..rhs{self[i]=T::zero();}}}implShrfor FormalPowerSerieswhere T:FormalPowerSeriesCoefficient{type Output=Self;fn shr(mut self,rhs:usize)->Self::Output{self.shr_assign(rhs);self}}implShlfor FormalPowerSerieswhere T:FormalPowerSeriesCoefficient{type Output=Self;fn shl(mut self,rhs:usize)->Self::Output{self.shl_assign(rhs);self}}implShrfor&FormalPowerSerieswhere T:FormalPowerSeriesCoefficient{type Output=FormalPowerSeries;fn shr(self,rhs:usize)->Self::Output{if self.length()<=rhs{Self::Output::zero()}else{let mut f=Self::Output::zeros(self.length()-rhs);for i in rhs..self.length(){f[i-rhs]=self[i].clone();}f}}}implShlfor&FormalPowerSerieswhere T:FormalPowerSeriesCoefficient{type Output=FormalPowerSeries;fn shl(self,rhs:usize)->Self::Output{let mut f=Self::Output::zeros(self.length()+rhs);for(i,x)in self.iter().cloned().enumerate().rev(){f[i+rhs]=x;}f}}}} // codesnip-guard: BarrettReduction #[cfg_attr(any(),rustfmt::skip)]pub use self::barrett_reduction::{BarrettReduction,Barrettable};#[cfg_attr(any(),rustfmt::skip)]mod barrett_reduction{use super::{One,Zero};use std::ops::{Add,Mul,Sub};#[derive(Debug,Clone,Copy)]pub struct BarrettReduction{m:T,im:T}implBarrettReductionwhere T:Barrettable{pub fn new(m:T)->Self{Self{m,im:T::inv_mod_approx(m)}}pub const fn new_with_im(m:T,im:T)->Self{Self{m,im}}pub const fn get_mod(&self)->T{self.m}pub fn div_rem(&self,a:T)->(T,T){T::barrett_reduce(a,self.m,self.im)}pub fn div(&self,a:T)->T{self.div_rem(a).0}pub fn rem(&self,a:T)->T{self.div_rem(a).1}}pub trait Barrettable:Sized+Copy+PartialOrd+Zero+One+Add+Sub+Mul{fn inv_mod_approx(m:Self)->Self;fn div_approx(self,im:Self)->Self;fn barrett_reduce(self,m:Self,im:Self)->(Self,Self){if m==Self::one(){return(self,Self::zero());}let q=self.div_approx(im);let r=self-q*m;if m<=r{(q+Self::one(),r-m)}else{(q,r)}}}impl Barrettable for u32{fn inv_mod_approx(m:Self)->Self{!0/m}fn div_approx(self,im:Self)->Self{((self as u64*im as u64)>>32)as u32}}impl Barrettable for u64{fn inv_mod_approx(m:Self)->Self{!0/m}fn div_approx(self,im:Self)->Self{((self as u128*im as u128)>>64)as u64}}impl Barrettable for u128{fn inv_mod_approx(m:Self)->Self{!0/m}fn div_approx(self,im:Self)->Self{const MASK64:u128=0xffff_ffff_ffff_ffff;let au=self>>64;let ad=self&MASK64;let imu=im>>64;let imd=im&MASK64;let res=au*imu;let x=((ad*imd)>>64)+au*imd;let(x,c)=x.overflowing_add(ad*imu);res+((c as u128)<<64)+(x>>64)}}} // codesnip-guard: Complex #[cfg_attr(any(),rustfmt::skip)]pub use self::complex::Complex;#[cfg_attr(any(),rustfmt::skip)]mod complex{use super::{Float,One,Scan,ScanSource,Zero};use std::{cmp::Ordering,iter::{Product,Sum},ops::{Add,AddAssign,Div,DivAssign,Mul,MulAssign,Neg,Sub,SubAssign}};#[derive(Clone,Copy,Debug,Default,PartialEq,Eq,PartialOrd,Ord,Hash)]pub struct Complex{pub re:T,pub im:T}implComplex{pub fn new(re:T,im:T)->Self{Self{re,im}}pub fn transpose(self)->Self{Self::new(self.im,self.re)}pub fn map(self,mut f:impl FnMut(T)->U)->Complex{Complex::new(f(self.re),f(self.im))}}implZero for Complexwhere T:Zero{fn zero()->Self{Self::new(T::zero(),T::zero())}}implOne for Complexwhere T:Zero+One{fn one()->Self{Self::new(T::one(),T::zero())}}implComplexwhere T:Zero+One{pub fn i()->Self{Self::new(T::zero(),T::one())}}implComplexwhere T:Neg{pub fn conjugate(self)->Self{Self::new(self.re,-self.im)}}implComplexwhere T:Mul{pub fn dot(self,rhs:Self)-><::Output as Add>::Output{self.re*rhs.re+self.im*rhs.im}}implComplexwhere T:Mul{pub fn cross(self,rhs:Self)-><::Output as Sub>::Output{self.re*rhs.im-self.im*rhs.re}}implComplexwhere T:Mul+Clone{pub fn norm(self)-><::Output as Add>::Output{self.re.clone()*self.re+self.im.clone()*self.im}}implComplexwhere T:Zero+Ord+Mul{pub fn cmp_by_arg(self,other:Self)->Ordering{fn pos(c:&Complex)->bool where T:Zero+Ord{let zero=T::zero();c.imComplexwhere T:Float{pub fn polar(r:T,theta:T)->Self{Self::new(r*theta.cos(),r*theta.sin())}pub fn primitive_nth_root_of_unity(n:T)->Self{let theta=T::TAU/n;Self::new(theta.cos(),theta.sin())}pub fn abs(self)->T{self.re.hypot(self.im)}pub fn unit(self)->Self{self/self.abs()}pub fn angle(self)->T{self.im.atan2(self.re)}}implAdd for Complexwhere T:Add{type Output=Complex<::Output>;fn add(self,rhs:Self)->Self::Output{Complex::new(self.re+rhs.re,self.im+rhs.im)}}implAddfor Complexwhere T:Add{type Output=Self;fn add(self,rhs:T)->Self::Output{Self::new(self.re+rhs,self.im)}}implSub for Complexwhere T:Sub{type Output=Complex<::Output>;fn sub(self,rhs:Self)->Self::Output{Complex::new(self.re-rhs.re,self.im-rhs.im)}}implSubfor Complexwhere T:Sub{type Output=Self;fn sub(self,rhs:T)->Self::Output{Self::new(self.re-rhs,self.im)}}implMul for Complexwhere T:Clone+Mul,::Output:Add+Sub{type Output=Complex;fn mul(self,rhs:Self)->Self::Output{Complex::new(self.re.clone()*rhs.re.clone()-self.im.clone()*rhs.im.clone(),self.re*rhs.im+self.im*rhs.re)}}implMulfor Complexwhere T:Clone+Mul{type Output=Complex<::Output>;fn mul(self,rhs:T)->Self::Output{Complex::new(self.re*rhs.clone(),self.im*rhs)}}implDiv for Complexwhere T:Clone+Add+Sub+Mul+Div{type Output=Complex<::Output>;fn div(self,rhs:Self)->Self::Output{let d=rhs.re.clone()*rhs.re.clone()+rhs.im.clone()*rhs.im.clone();Complex::new((self.re.clone()*rhs.re.clone()+self.im.clone()*rhs.im.clone())/d.clone(),(self.im*rhs.re-self.re*rhs.im)/d)}}implDivfor Complexwhere T:Clone+Div{type Output=Complex<::Output>;fn div(self,rhs:T)->Self::Output{Complex::new(self.re/rhs.clone(),self.im/rhs)}}implNeg for Complexwhere T:Neg{type Output=Complex<::Output>;fn neg(self)->Self::Output{Complex::new(-self.re,-self.im)}}macro_rules!impl_complex_ref_binop{(impl<$T:ident>$imp:ident$method:ident($l:ty,$r:ty)where$($w:ident)*)=>{impl<$T>$imp<$r>for&$l where$T:Clone$(+$w)*,{type Output=<$l as$imp<$r>>::Output;fn$method(self,rhs:$r)-><$l as$imp<$r>>::Output{$imp::$method(self.clone(),rhs)}}impl<$T>$imp<&$r>for$l where$T:Clone$(+$w)*,{type Output=<$l as$imp<$r>>::Output;fn$method(self,rhs:&$r)-><$l as$imp<$r>>::Output{$imp::$method(self,rhs.clone())}}impl<$T>$imp<&$r>for&$l where$T:Clone$(+$w)*,{type Output=<$l as$imp<$r>>::Output;fn$method(self,rhs:&$r)-><$l as$imp<$r>>::Output{$imp::$method(self.clone(),rhs.clone())}}};}impl_complex_ref_binop!(implAdd add(Complex,Complex)where Add);impl_complex_ref_binop!(implAdd add(Complex,T)where Add);impl_complex_ref_binop!(implSub sub(Complex,Complex)where Sub);impl_complex_ref_binop!(implSub sub(Complex,T)where Sub);impl_complex_ref_binop!(implMul mul(Complex,Complex)where Add Sub Mul);impl_complex_ref_binop!(implMul mul(Complex,T)where Mul);impl_complex_ref_binop!(implDiv div(Complex,Complex)where Add Sub Mul Div);impl_complex_ref_binop!(implDiv div(Complex,T)where Div);macro_rules!impl_complex_ref_unop{(impl<$T:ident>$imp:ident$method:ident($t:ty)where$($w:ident)*)=>{impl<$T>$imp for&$t where$T:Clone$(+$w)*,{type Output=<$t as$imp>::Output;fn$method(self)-><$t as$imp>::Output{$imp::$method(self.clone())}}};}impl_complex_ref_unop!(implNeg neg(Complex)where Neg);macro_rules!impl_complex_op_assign{(impl<$T:ident>$imp:ident$method:ident($l:ty,$r:ty)$fromimp:ident$frommethod:ident where$($w:ident)*)=>{impl<$T>$imp<$r>for$l where$T:Clone$(+$w)*,{fn$method(&mut self,rhs:$r){*self=$fromimp::$frommethod(self.clone(),rhs);}}impl<$T>$imp<&$r>for$l where$T:Clone$(+$w)*,{fn$method(&mut self,rhs:&$r){$imp::$method(self,rhs.clone());}}};}impl_complex_op_assign!(implAddAssign add_assign(Complex,Complex)Add add where Add);impl_complex_op_assign!(implAddAssign add_assign(Complex,T)Add add where Add);impl_complex_op_assign!(implSubAssign sub_assign(Complex,Complex)Sub sub where Sub);impl_complex_op_assign!(implSubAssign sub_assign(Complex,T)Sub sub where Sub);impl_complex_op_assign!(implMulAssign mul_assign(Complex,Complex)Mul mul where Add Sub Mul);impl_complex_op_assign!(implMulAssign mul_assign(Complex,T)Mul mul where Mul);impl_complex_op_assign!(implDivAssign div_assign(Complex,Complex)Div div where Add Sub Mul Div);impl_complex_op_assign!(implDivAssign div_assign(Complex,T)Div div where Div);macro_rules!impl_complex_fold{(impl<$T:ident>$imp:ident$method:ident($t:ty)$identimp:ident$identmethod:ident$fromimp:ident$frommethod:ident where$($w:ident)*$(+$x:ident)*)=>{impl<$T>$imp for$t where$T:$identimp$(+$w)*$(+$x)*,{fn$method>(iter:I)->Self{iter.fold(::$identmethod(),$fromimp::$frommethod)}}impl<'a,$T:'a>$imp<&'a$t>for$t where$T:Clone+$identimp$(+$w)*$(+$x)*,{fn$method>(iter:I)->Self{iter.fold(::$identmethod(),$fromimp::$frommethod)}}};}impl_complex_fold!(implSum sum(Complex)Zero zero Add add where Add);impl_complex_fold!(implProduct product(Complex)One one Mul mul where Add Sub Mul+Zero+Clone);implScan for Complex{type Output=Complex<::Output>;fn scan(iter:&mut I)->Option{Some(Complex::new(::scan(iter)?,::scan(iter)?))}}} // codesnip-guard: ConvolveSteps #[cfg_attr(any(),rustfmt::skip)]pub use self::convolve_steps::ConvolveSteps;#[cfg_attr(any(),rustfmt::skip)]mod convolve_steps{pub trait ConvolveSteps{#[doc=" Whether transform multiplication computes modulo x^n - 1 in the coefficient ring."]const CYCLIC:bool=false;type T;type F;fn length(t:&Self::T)->usize;fn transform(t:Self::T,len:usize)->Self::F;fn inverse_transform(f:Self::F,len:usize)->Self::T;fn multiply(f:&mut Self::F,g:&Self::F);fn square(t:Self::T,len:usize)->Self::T where Self::T:Clone{let mut f=Self::transform(t.clone(),len);let g=Self::transform(t,len);Self::multiply(&mut f,&g);Self::inverse_transform(f,len)}fn convolve(a:Self::T,b:Self::T)->Self::T{let len=(Self::length(&a)+Self::length(&b)).saturating_sub(1);let mut a=Self::transform(a,len);let b=Self::transform(b,len);Self::multiply(&mut a,&b);Self::inverse_transform(a,len)}}} // codesnip-guard: MInt #[cfg_attr(any(),rustfmt::skip)]pub mod mint_basic{use super::*;use std::{cell::UnsafeCell,mem::swap};#[macro_export]macro_rules!define_basic_mintbase{($name:ident,$m:expr,$basety:ty,$signedty:ty,$upperty:ty,[$($unsigned:ty),*],[$($signed:ty),*])=>{pub enum$name{}impl MIntBase for$name{type Inner=$basety;#[inline]fn get_mod()->Self::Inner{$m}#[inline]fn mod_zero()->Self::Inner{0}#[inline]fn mod_one()->Self::Inner{(Self::get_mod()!=1)as$basety}#[inline]fn mod_add(x:Self::Inner,y:Self::Inner)->Self::Inner{let m=Self::get_mod();let(z,borrow)=x.overflowing_sub(m-y);if borrow{z.wrapping_add(m)}else{z}}#[inline]fn mod_sub(x:Self::Inner,y:Self::Inner)->Self::Inner{if xSelf::Inner{$name::rem(x as$upperty*y as$upperty)as$basety}#[inline]fn mod_div(x:Self::Inner,y:Self::Inner)->Self::Inner{Self::mod_mul(x,Self::mod_inv(y))}#[inline]fn mod_neg(x:Self::Inner)->Self::Inner{if x==0{0}else{Self::get_mod()-x}}fn mod_inv(x:Self::Inner)->Self::Inner{let(mut a,mut b)=(x,Self::get_mod());let(mut u,mut v)=(1,0);let mut negative=true;while a!=0{let k=b/a;v+=k*u;b-=k*a;swap(&mut u,&mut v);swap(&mut b,&mut a);negative=!negative;}if negative{Self::mod_neg(v)}else{v}}}$(impl MIntConvert<$unsigned>for$name{#[inline]fn from(x:$unsigned)->Self::Inner{(x%::get_mod()as$unsigned)as$basety}#[inline]fn into(x:Self::Inner)->$unsigned{x as$unsigned}#[inline]fn mod_into()->$unsigned{::get_mod()as$unsigned}})*$(impl MIntConvert<$signed>for$name{#[inline]fn from(x:$signed)->Self::Inner{let modulus=(::get_mod()as$signed).cast_unsigned();let value=(x.unsigned_abs()%modulus)as$basety;if x<0{::mod_neg(value)}else{value}}#[inline]fn into(x:Self::Inner)->$signed{x as$signed}#[inline]fn mod_into()->$signed{::get_mod()as$signed}})*};}#[macro_export]macro_rules!define_basic_mint32{($([$name:ident,$m:expr,$mint_name:ident]),*)=>{$(define_basic_mintbase!($name,$m,u32,i32,u64,[u32,u64,u128,usize],[i32,i64,i128,isize]);impl$name{fn rem(x:u64)->u64{x%$m}}pub type$mint_name=MInt<$name>;)*};}thread_local!(static DYN_MODULUS_U32:UnsafeCell> =const{UnsafeCell::new(BarrettReduction::::new_with_im(1_000_000_007,!0/1_000_000_007))});impl DynModuloU32{pub fn set_mod(m:u32){DYN_MODULUS_U32.with(|cell|unsafe{*cell.get()=BarrettReduction::::new(m as u64)});}fn rem(x:u64)->u64{DYN_MODULUS_U32.with(|cell|unsafe{(*cell.get()).rem(x)})}}impl DynMIntU32{pub fn set_mod(m:u32){DynModuloU32::set_mod(m)}}thread_local!(static DYN_MODULUS_U64:UnsafeCell> =const{UnsafeCell::new(BarrettReduction::::new_with_im(1_000_000_007,!0/1_000_000_007))});impl DynModuloU64{pub fn set_mod(m:u64){DYN_MODULUS_U64.with(|cell|unsafe{*cell.get()=BarrettReduction::::new(m as u128)})}fn rem(x:u128)->u128{DYN_MODULUS_U64.with(|cell|unsafe{(*cell.get()).rem(x)})}}impl DynMIntU64{pub fn set_mod(m:u64){DynModuloU64::set_mod(m)}}define_basic_mint32!([Modulo998244353,998_244_353,MInt998244353],[Modulo1000000007,1_000_000_007,MInt1000000007],[Modulo1000000009,1_000_000_009,MInt1000000009]);define_basic_mintbase!(DynModuloU32,DYN_MODULUS_U32.with(|cell|unsafe{(*cell.get()).get_mod()as u32}),u32,i32,u64,[u32,u64,u128,usize],[i32,i64,i128,isize]);pub type DynMIntU32=MInt;define_basic_mintbase!(DynModuloU64,DYN_MODULUS_U64.with(|cell|unsafe{(*cell.get()).get_mod()as u64}),u64,i64,u128,[u64,u128,usize],[i64,i128,isize]);pub type DynMIntU64=MInt;pub struct Modulo2;impl MIntBase for Modulo2{type Inner=u32;#[inline]fn get_mod()->Self::Inner{2}#[inline]fn mod_zero()->Self::Inner{0}#[inline]fn mod_one()->Self::Inner{1}#[inline]fn mod_add(x:Self::Inner,y:Self::Inner)->Self::Inner{x^y}#[inline]fn mod_sub(x:Self::Inner,y:Self::Inner)->Self::Inner{x^y}#[inline]fn mod_mul(x:Self::Inner,y:Self::Inner)->Self::Inner{x&y}#[inline]fn mod_div(x:Self::Inner,y:Self::Inner)->Self::Inner{assert_ne!(y,0);x}#[inline]fn mod_neg(x:Self::Inner)->Self::Inner{x}#[inline]fn mod_inv(x:Self::Inner)->Self::Inner{assert_ne!(x,0);x}#[inline]fn mod_pow(x:Self::Inner,y:usize)->Self::Inner{if y==0{1}else{x}}}macro_rules!impl_to_mint_base_for_modulo2{($name:ident,$basety:ty,[$($t:ty),*])=>{$(impl MIntConvert<$t>for$name{#[inline]fn from(x:$t)->Self::Inner{(x&1)as$basety}#[inline]fn into(x:Self::Inner)->$t{x as$t}#[inline]fn mod_into()->$t{2}})*};}impl_to_mint_base_for_modulo2!(Modulo2,u32,[u8,u16,u32,u64,u128,usize,i8,i16,i32,i64,i128,isize]);pub type MInt2=MInt;} // codesnip-guard: MIntBase #[cfg_attr(any(),rustfmt::skip)]pub use mint_base::{MInt,MIntBase,MIntConvert};#[cfg_attr(any(),rustfmt::skip)]mod mint_base{use super::*;use std::{fmt::{self,Debug,Display},hash::{Hash,Hasher},iter::{Product,Sum},marker::PhantomData,ops::{Add,AddAssign,Div,DivAssign,Mul,MulAssign,Neg,Sub,SubAssign},str::FromStr};#[repr(transparent)]pub struct MIntwhere M:MIntBase{x:M::Inner,_marker:PhantomDataM>}pub trait MIntConvert::Inner>:MIntBase{fn from(x:T)->::Inner;fn into(x:::Inner)->T;fn mod_into()->T;}pub trait MIntBase{type Inner:Sized+Copy+Eq+Debug+Hash;fn get_mod()->Self::Inner;fn mod_zero()->Self::Inner;fn mod_one()->Self::Inner;fn mod_add(x:Self::Inner,y:Self::Inner)->Self::Inner;fn mod_sub(x:Self::Inner,y:Self::Inner)->Self::Inner;fn mod_mul(x:Self::Inner,y:Self::Inner)->Self::Inner;fn mod_div(x:Self::Inner,y:Self::Inner)->Self::Inner;fn mod_neg(x:Self::Inner)->Self::Inner;fn mod_inv(x:Self::Inner)->Self::Inner;fn mod_pow(x:Self::Inner,y:usize)->Self::Inner{let(mut x,mut y,mut z)=(x,y,Self::mod_one());while y>0{if y&1==1{z=Self::mod_mul(z,x);}x=Self::mod_mul(x,x);y>>=1;}z}fn mod_inner(x:Self::Inner)->Self::Inner{x}}implMIntwhere M:MIntConvert{#[inline]pub fn new(x:M::Inner)->Self{Self::new_unchecked(>::from(x))}}implMIntwhere M:MIntBase{#[inline]pub const fn new_unchecked(x:M::Inner)->Self{Self{x,_marker:PhantomData}}#[inline]pub fn get_mod()->M::Inner{M::get_mod()}#[inline]pub fn pow(self,y:usize)->Self{Self::new_unchecked(M::mod_pow(self.x,y))}#[inline]pub fn inv(self)->Self{Self::new_unchecked(M::mod_inv(self.x))}#[inline]pub fn inner(self)->M::Inner{M::mod_inner(self.x)}}implClone for MIntwhere M:MIntBase{#[inline]fn clone(&self)->Self{*self}}implCopy for MIntwhere M:MIntBase{}implDebug for MIntwhere M:MIntBase{fn fmt(&self,f:&mut fmt::Formatter<'_>)->fmt::Result{Debug::fmt(&self.inner(),f)}}implDefault for MIntwhere M:MIntBase{#[inline]fn default()->Self{::zero()}}implPartialEq for MIntwhere M:MIntBase{#[inline]fn eq(&self,other:&Self)->bool{PartialEq::eq(&self.x,&other.x)}}implEq for MIntwhere M:MIntBase{}implHash for MIntwhere M:MIntBase{#[inline]fn hash(&self,state:&mut H){Hash::hash(&self.x,state)}}macro_rules!impl_mint_from{($($t:ty),*)=>{$(implFrom<$t>for MIntwhere M:MIntConvert<$t>,{#[inline]fn from(x:$t)->Self{Self::new_unchecked(>::from(x))}}implFrom>for$t where M:MIntConvert<$t>,{#[inline]fn from(x:MInt)->$t{>::into(x.x)}})*};}impl_mint_from!(u8,u16,u32,u64,u128,usize,i8,i16,i32,i64,i128,isize);implZero for MIntwhere M:MIntBase{#[inline]fn zero()->Self{Self::new_unchecked(M::mod_zero())}}implOne for MIntwhere M:MIntBase{#[inline]fn one()->Self{Self::new_unchecked(M::mod_one())}}implAdd for MIntwhere M:MIntBase{type Output=Self;#[inline]fn add(self,rhs:Self)->Self::Output{Self::new_unchecked(M::mod_add(self.x,rhs.x))}}implSub for MIntwhere M:MIntBase{type Output=Self;#[inline]fn sub(self,rhs:Self)->Self::Output{Self::new_unchecked(M::mod_sub(self.x,rhs.x))}}implMul for MIntwhere M:MIntBase{type Output=Self;#[inline]fn mul(self,rhs:Self)->Self::Output{Self::new_unchecked(M::mod_mul(self.x,rhs.x))}}implDiv for MIntwhere M:MIntBase{type Output=Self;#[inline]fn div(self,rhs:Self)->Self::Output{Self::new_unchecked(M::mod_div(self.x,rhs.x))}}implNeg for MIntwhere M:MIntBase{type Output=Self;#[inline]fn neg(self)->Self::Output{Self::new_unchecked(M::mod_neg(self.x))}}implSum for MIntwhere M:MIntBase{#[inline]fn sum>(iter:I)->Self{iter.fold(::zero(),Add::add)}}implProduct for MIntwhere M:MIntBase{#[inline]fn product>(iter:I)->Self{iter.fold(::one(),Mul::mul)}}impl<'a,M:'a>Sum<&'a MInt>for MIntwhere M:MIntBase{#[inline]fn sum>(iter:I)->Self{iter.fold(::zero(),Add::add)}}impl<'a,M:'a>Product<&'a MInt>for MIntwhere M:MIntBase{#[inline]fn product>(iter:I)->Self{iter.fold(::one(),Mul::mul)}}implDisplay for MIntwhere M:MIntBase{fn fmt(&self,f:&mut fmt::Formatter<'_>)->Result<(),fmt::Error>{write!(f,"{}",self.inner())}}implFromStr for MIntwhere M:MIntConvert+MIntBase{type Err=::Err;#[inline]fn from_str(s:&str)->Result{s.parse::().map(Self::new)}}implScan for MIntwhere M:MIntConvert+MIntBase>{type Output=Self;#[inline]fn scan(source:&mut S)->Option{M::Inner::scan(source).map(Self::new)}fn scan_vec(source:&mut S,size:usize,capacity:usize)->Vec{let mut values:Vec<_>=source.scan_vec::(size,capacity).into_iter().map(Self::new).collect();values.reserve_exact(capacity.saturating_sub(values.len()));values}}macro_rules!impl_mint_ref_binop{($imp:ident,$method:ident,$t:ty)=>{impl$imp<$t>for&$t where M:MIntBase,{type Output=<$t as$imp<$t>>::Output;#[inline]fn$method(self,other:$t)-><$t as$imp<$t>>::Output{$imp::$method(*self,other)}}impl$imp<&$t>for$t where M:MIntBase,{type Output=<$t as$imp<$t>>::Output;#[inline]fn$method(self,other:&$t)-><$t as$imp<$t>>::Output{$imp::$method(self,*other)}}impl$imp<&$t>for&$t where M:MIntBase,{type Output=<$t as$imp<$t>>::Output;#[inline]fn$method(self,other:&$t)-><$t as$imp<$t>>::Output{$imp::$method(*self,*other)}}};}impl_mint_ref_binop!(Add,add,MInt);impl_mint_ref_binop!(Sub,sub,MInt);impl_mint_ref_binop!(Mul,mul,MInt);impl_mint_ref_binop!(Div,div,MInt);macro_rules!impl_mint_ref_unop{($imp:ident,$method:ident,$t:ty)=>{impl$imp for&$t where M:MIntBase,{type Output=<$t as$imp>::Output;#[inline]fn$method(self)-><$t as$imp>::Output{$imp::$method(*self)}}};}impl_mint_ref_unop!(Neg,neg,MInt);macro_rules!impl_mint_ref_op_assign{($imp:ident,$method:ident,$t:ty,$fromimp:ident,$frommethod:ident)=>{impl$imp<$t>for$t where M:MIntBase,{#[inline]fn$method(&mut self,rhs:$t){*self=$fromimp::$frommethod(*self,rhs);}}impl$imp<&$t>for$t where M:MIntBase,{#[inline]fn$method(&mut self,other:&$t){$imp::$method(self,*other);}}};}impl_mint_ref_op_assign!(AddAssign,add_assign,MInt,Add,add);impl_mint_ref_op_assign!(SubAssign,sub_assign,MInt,Sub,sub);impl_mint_ref_op_assign!(MulAssign,mul_assign,MInt,Mul,mul);impl_mint_ref_op_assign!(DivAssign,div_assign,MInt,Div,div);} // codesnip-guard: NumberTheoreticTransform #[cfg_attr(any(),rustfmt::skip)]pub use self::number_theoretic_transform::{Convolve,Convolve998244353,MIntConvolve,NttReuse,U64Convolve};#[cfg_attr(any(),rustfmt::skip)]#[cfg(target_arch="x86_64")]mod mint_fft_convolve{#![allow(unsafe_op_in_unsafe_fn)]#[cfg(target_arch="x86_64")]use super::{MInt,MIntConvert,advise_huge_pages,partial_fft::{avx2::{Complex4,complex_mul,fft_soa,ifft_soa,load,multiply_accumulate,store},eval_twiddle,with_roots}};#[cfg(target_arch="x86_64")]use std::arch::x86_64::*;#[target_feature(enable="avx2,fma")]#[inline]#[cfg(target_arch="x86_64")]unsafe fn round4(value:&[f64;4])->[i64;4]{let magic=_mm256_set1_pd((3i64<<51)as f64);let rounded=_mm256_sub_epi64(_mm256_castpd_si256(_mm256_add_pd(_mm256_load_pd(value.as_ptr()),magic)),_mm256_castpd_si256(magic));let mut result=[0;4];_mm256_storeu_si256(result.as_mut_ptr().cast(),rounded);result}#[target_feature(enable="avx2,fma")]#[inline]#[cfg(target_arch="x86_64")]unsafe fn reduce_mod4(value:__m256d,modulus:__m256d,inverse:__m256d)->__m256d{let quotient=_mm256_round_pd::<{_MM_FROUND_TO_NEAREST_INT|_MM_FROUND_NO_EXC}>(_mm256_mul_pd(value,inverse));_mm256_fnmadd_pd(quotient,modulus,value)}#[target_feature(enable="avx2,fma")]#[cfg(target_arch="x86_64")]unsafe fn split_coefficients(values:Vec>,n:usize,modulus:i64,split:i64)->(Vec,Vec)where M:MIntConvert+MIntConvert{let mut low=Vec::::with_capacity(n/4);let mut high=Vec::::with_capacity(n/4);advise_huge_pages(&mut low);advise_huge_pages(&mut high);let divisor=_mm256_set1_pd(split as f64);let split=_mm256_set1_pd(split as f64);for i in(0..values.len()).step_by(4){let mut centered=[0i32;4];for lane in 0..4 .min(values.len()-i){let mut value=u32::from(values[i+lane])as i64;if value*2>modulus{value-=modulus;}centered[lane]=value as i32;}let value=_mm256_cvtepi32_pd(_mm_loadu_si128(centered.as_ptr().cast()));let upper=_mm256_round_pd::<{_MM_FROUND_TO_NEAREST_INT|_MM_FROUND_NO_EXC}>(_mm256_div_pd(value,divisor));let lower=_mm256_fnmadd_pd(upper,split,value);if i>2].im.as_mut_ptr(),lower);_mm256_store_pd(high[(i-n)>>2].im.as_mut_ptr(),upper);}}low.resize(n/4,Complex4::default());high.resize(n/4,Complex4::default());(low,high)}#[target_feature(enable="avx2,fma")]#[cfg(target_arch="x86_64")]unsafe fn dot_soa(a0:&mut[Complex4],a1:&mut[Complex4],b0:&mut[Complex4],b1:&[Complex4]){let n=a0.len()*4;with_roots(n/2,|cache|{for i in 0..a0.len(){let(mut cr,mut ci)=load(&b0[i]);let(mut dr,mut di)=load(&b1[i]);let mut c0r=_mm256_setzero_pd();let mut c0i=_mm256_setzero_pd();let mut c1r=_mm256_setzero_pd();let mut c1i=_mm256_setzero_pd();let mut c2r=_mm256_setzero_pd();let mut c2i=_mm256_setzero_pd();let w=eval_twiddle(cache,1,a0.len(),i);let wr=_mm256_setr_pd(w.re,1.0,1.0,1.0);let wi=_mm256_setr_pd(w.im,0.0,0.0,0.0);for lane in 0..4{let ar=_mm256_set1_pd(a0[i].re[lane]);let ai=_mm256_set1_pd(a0[i].im[lane]);let br=_mm256_set1_pd(a1[i].re[lane]);let bi=_mm256_set1_pd(a1[i].im[lane]);multiply_accumulate(&mut c0r,&mut c0i,ar,ai,cr,ci);multiply_accumulate(&mut c1r,&mut c1i,ar,ai,dr,di);multiply_accumulate(&mut c1r,&mut c1i,br,bi,cr,ci);multiply_accumulate(&mut c2r,&mut c2i,br,bi,dr,di);if lane!=3{cr=_mm256_permute4x64_pd::<0x93>(cr);ci=_mm256_permute4x64_pd::<0x93>(ci);dr=_mm256_permute4x64_pd::<0x93>(dr);di=_mm256_permute4x64_pd::<0x93>(di);(cr,ci)=complex_mul(cr,ci,wr,wi);(dr,di)=complex_mul(dr,di,wr,wi);}}store(&mut a0[i],c0r,c0i);store(&mut a1[i],c1r,c1i);store(&mut b0[i],c2r,c2i);}});}#[target_feature(enable="avx2,fma")]#[cfg(target_arch="x86_64")]unsafe fn split_u64_coefficients(values:&[u64],n:usize)->[Vec;5]{let mut result:[Vec;5]=std::array::from_fn(|_|{let mut part=Vec::with_capacity(n/4);advise_huge_pages(&mut part);part});for(i,chunk)in values.chunks(4).enumerate(){let mut parts=[Complex4::default();5];for(lane,mut value)in chunk.iter().copied().enumerate(){for part in&mut parts{let digit=((value<<51)as i64)>>51;value=(value>>13).wrapping_add(u64::from(digit<0));part.re[lane]=digit as f64;}}for(result,part)in result.iter_mut().zip(parts){if i(a:&mut[Vec;5],b:&[Vec;5]){let n=a[0].len()*4;with_roots(n/2,|cache|{for block in 0..a[0].len(){let mut br=[_mm256_setzero_pd();PARTS];let mut bi=br;let mut rr=[_mm256_setzero_pd();5];let mut ri=rr;for part in 0..PARTS{(br[part],bi[part])=load(&b[part][block]);}let w=eval_twiddle(cache,1,a[0].len(),block);let wr=_mm256_setr_pd(w.re,1.0,1.0,1.0);let wi=_mm256_setr_pd(w.im,0.0,0.0,0.0);macro_rules!accumulate_lane{($lane:expr)=>{{for right in(0..PARTS).rev(){for left in 0..PARTS.min(5-right){let ar=_mm256_set1_pd(a[left][block].re[$lane]);let ai=_mm256_set1_pd(a[left][block].im[$lane]);multiply_accumulate(&mut rr[left+right],&mut ri[left+right],ar,ai,br[right],bi[right],);}if$lane!=3{br[right]=_mm256_permute4x64_pd::<0x93>(br[right]);bi[right]=_mm256_permute4x64_pd::<0x93>(bi[right]);(br[right],bi[right])=complex_mul(br[right],bi[right],wr,wi);}}}};}accumulate_lane!(0);accumulate_lane!(1);accumulate_lane!(2);accumulate_lane!(3);for part in 0..5{store(&mut a[part][block],rr[part],ri[part]);}}});}#[target_feature(enable="avx2,fma")]#[cfg(target_arch="x86_64")]pub unsafe fn convolve_u64_avx2(a:Vec,b:Vec)->Vec{let len=a.len()+b.len()-1;let n=len.next_power_of_two()/2;let parts=|values:&[u64]|{let maximum=values.iter().try_fold(0,|maximum,&value|{(value<=4095*(1+(1u64<<13)+(1u64<<26)+(1u64<<39))).then_some(maximum.max(value))}).unwrap_or(u64::MAX);let mut limit=(1u64<<12)-1;let mut parts=1;while parts<5&&maximum>limit{limit=(limit<<13)+((1<<12)-1);parts+=1;}parts};let a_parts=parts(&a);let mut fa=split_u64_coefficients(&a,n);drop(a);let b_parts=parts(&b);let mut fb=split_u64_coefficients(&b,n);drop(b);let shared=a_parts.min(b_parts).min(3);for part in 0..shared{fft_soa(&mut fa[part]);fft_soa(&mut fb[part]);}for part in&mut fa[shared..a_parts]{fft_soa(part);}for part in&mut fb[shared..b_parts]{fft_soa(part);}match a_parts.max(b_parts){1=>dot_u64_soa::<1>(&mut fa,&fb),2=>dot_u64_soa::<2>(&mut fa,&fb),3=>dot_u64_soa::<3>(&mut fa,&fb),4=>dot_u64_soa::<4>(&mut fa,&fb),_=>dot_u64_soa::<5>(&mut fa,&fb),}drop(fb);for part in&mut fa[..5 .min(a_parts+b_parts-1)]{ifft_soa(part);}let mut result=vec![0;len];for(block,_)in fa[0].iter().enumerate(){let real:[[i64;4];5]=std::array::from_fn(|part|round4(&fa[part][block].re));let imag:[[i64;4];5]=std::array::from_fn(|part|round4(&fa[part][block].im));for lane in 0..4{let i=block*4+lane;if i(a:Vec>,b:Vec>)->Vec>where M:MIntConvert+MIntConvert{let len=a.len()+b.len()-1;let n=len.next_power_of_two()/2;let modulus=>::mod_into()as i64;let split=(modulus as f64).sqrt()as i64+1;let(mut a0,mut a1)=split_coefficients(a,n,modulus,split);let(mut b0,mut b1)=split_coefficients(b,n,modulus,split);fft_soa(&mut a0);fft_soa(&mut a1);fft_soa(&mut b0);fft_soa(&mut b1);dot_soa(&mut a0,&mut a1,&mut b0,&b1);drop(b1);ifft_soa(&mut a0);ifft_soa(&mut a1);ifft_soa(&mut b0);let split2=(split*split%modulus)as f64;let split=_mm256_set1_pd(split as f64);let split2=_mm256_set1_pd(split2);let inverse=_mm256_set1_pd(1.0/modulus as f64);let modulus=_mm256_set1_pd(modulus as f64);let magic=_mm256_set1_pd((3i64<<51)as f64);let mut result=vec![MInt::::from(0u32);len];for(block,((a0,a1),b0))in a0.iter().zip(&a1).zip(&b0).enumerate(){for(part,(a0,a1,b0))in[(&a0.re,&a1.re,&b0.re),(&a0.im,&a1.im,&b0.im)].into_iter().enumerate(){let a0=_mm256_round_pd::<{_MM_FROUND_TO_NEAREST_INT|_MM_FROUND_NO_EXC}>(_mm256_load_pd(a0.as_ptr()));let a1=_mm256_round_pd::<{_MM_FROUND_TO_NEAREST_INT|_MM_FROUND_NO_EXC}>(_mm256_load_pd(a1.as_ptr()));let b0=_mm256_round_pd::<{_MM_FROUND_TO_NEAREST_INT|_MM_FROUND_NO_EXC}>(_mm256_load_pd(b0.as_ptr()));let a0=reduce_mod4(a0,modulus,inverse);let a1=reduce_mod4(a1,modulus,inverse);let b0=reduce_mod4(b0,modulus,inverse);let value=_mm256_fmadd_pd(b0,split2,_mm256_fmadd_pd(a1,split,a0));let value=reduce_mod4(value,modulus,inverse);let value=_mm256_add_pd(value,_mm256_and_pd(_mm256_cmp_pd::<_CMP_LT_OQ>(value,_mm256_setzero_pd()),modulus));let value=_mm256_sub_epi64(_mm256_castpd_si256(_mm256_add_pd(value,magic)),_mm256_castpd_si256(magic));let mut lanes=[0i64;4];_mm256_storeu_si256(lanes.as_mut_ptr().cast(),value);for(lane,value)in lanes.into_iter().enumerate(){let i=block*4+lane+part*n;if i>::mod_into();result[i]=MInt::::from(if valueSimdBackend{if width<8||!width.is_power_of_two()||len(PhantomDataM>);pub type Convolve998244353=Convolve;#[doc=" Raw transforms require each integer coefficient reconstructed by CRT to be below"]#[doc=" the product of the three NTT moduli. `convolve` splits products exceeding this bound."]pub type MIntConvolve=Convolve<(M,(Modulo167772161,Modulo469762049,Modulo754974721))>;#[doc=" Convolution modulo 2^64. Multiply only freshly transformed operands; reconstruct"]#[doc=" and transform again before multiplying another factor."]pub type U64Convolve=Convolve<(u64,(Modulo167772161,Modulo469762049,Modulo754974721))>;impl Montgomery32NttModulus for Modulo167772161{}impl Montgomery32NttModulus for Modulo469762049{}impl Montgomery32NttModulus for Modulo754974721{}impl Montgomery32NttModulus for Modulo998244353{}const fn reduce(z:u64,p:u32,r:u32)->u32{let mut z=((z+r.wrapping_mul(z as u32)as u64*p as u64)>>32)as u32;if z>=p{z-=p;}z}const fn mod_mul(x:u32,y:u32,p:u32,r:u32)->u32{reduce(x as u64*y as u64,p,r)}const fn mod_pow(mut x:u32,mut y:u32,p:u32,r:u32,mut z:u32)->u32{while y>0{if y&1==1{z=mod_mul(z,x,p,r);}x=mod_mul(x,x,p,r);y>>=1;}z}pub trait Montgomery32NttModulus:Sized+MontgomeryReduction32{const PRIMITIVE_ROOT:u32={let mut g=3u32;loop{let mut ok=true;let mut d=1u32;while d*d();}#[derive(Debug,PartialEq)]pub struct NttInfo{root:[u32;32],inv_root:[u32;32],rate3_packed:[[u32;8];32],inv_rate3_packed:[[u32;8];32]}impl NttInfo{const fn new()->Self where M:Montgomery32NttModulus{let mut root=[0;32];let mut inv_root=[0;32];let mut rate3_packed=[[0;8];32];let mut inv_rate3_packed=[[0;8];32];let rank=M::RANK as usize;let g=reduce(M::PRIMITIVE_ROOT as u64*M::N2 as u64,M::MOD,M::R);root[rank]=mod_pow(g,(M::MOD-1)>>rank,M::MOD,M::R,M::N1);inv_root[rank]=mod_pow(root[rank],M::MOD-2,M::MOD,M::R,M::N1);let mut i=rank-1;loop{root[i]=mod_mul(root[i+1],root[i+1],M::MOD,M::R);inv_root[i]=mod_mul(inv_root[i+1],inv_root[i+1],M::MOD,M::R);if i==0{break;}i-=1;}let(mut i,mut prod,mut inv_prod)=(0,M::N1,M::N1);while i(x:u32,y:u32)->u32 where M:Montgomery32NttModulus{let modulus=if M::MOD=modulus{sum-modulus}else{sum}}#[inline]fn sub_scalar(x:u32,y:u32)->u32 where M:Montgomery32NttModulus{let modulus=if M::MOD(x:u32,y:u32)->u32 where M:Montgomery32NttModulus{if M::MOD>32)as u32}else{M::mod_mul(x,y)}}fn ntt_scalar(a:&mut[MInt])where M:Montgomery32NttModulus{ntt_batch_scalar(a,1);}fn ntt_batch(a:&mut[MInt],width:usize)where M:Montgomery32NttModulus{#[cfg(target_arch="x86_64")]{match batch_ntt_simd_backend(a.len(),width){SimdBackend::Avx512=>{unsafe{ntt_simd::ntt_batch_avx512(a,width)};return;}SimdBackend::Avx2=>{unsafe{ntt_simd::ntt_batch_avx2::<_,false>(a,width)};return;}SimdBackend::Scalar=>{}}}ntt_batch_scalar(a,width);}fn ntt_batch_scalar(a:&mut[MInt],width:usize)where M:Montgomery32NttModulus{let n=a.len()/width;if n<=1{return;}let mut v=n/2;if n.trailing_zeros()&1==1{let(l,r)=a.split_at_mut(v*width);for(x0,x1)in l.iter_mut().zip(r){let a0=*x0;let a1=*x1;*x0=a0+a1;*x1=a0-a1;}v>>=1;}let imag=MInt::::new_unchecked(M::INFO.root[2]);while v>1{let mut w1=MInt::::one();let mut w2=w1;let mut w3=w1;for(s,a)in a.chunks_exact_mut((v<<1)*width).enumerate(){let(l,r)=a.split_at_mut(v*width);let(ll,lr)=l.split_at_mut((v>>1)*width);let(rl,rr)=r.split_at_mut((v>>1)*width);for(((x0,x1),x2),x3)in ll.iter_mut().zip(lr).zip(rl).zip(rr){let a0=*x0;let a1=*x1*w1;let a2=*x2*w2;let a3=*x3*w3;let a0pa2=a0+a2;let a0na2=a0-a2;let a1pa3=a1+a3;let a1na3imag=(a1-a3)*imag;*x0=a0pa2+a1pa3;*x1=a0pa2-a1pa3;*x2=a0na2+a1na3imag;*x3=a0na2-a1na3imag;}let rate=&M::INFO.rate3_packed[s.trailing_ones()as usize];w1*=MInt::::new_unchecked(rate[1]);w2*=MInt::::new_unchecked(rate[3]);w3*=MInt::::new_unchecked(rate[5]);}v>>=2;}}fn intt_scalar(a:&mut[MInt])where M:Montgomery32NttModulus{intt_batch_scalar(a,1);}fn intt_batch(a:&mut[MInt],width:usize)where M:Montgomery32NttModulus{#[cfg(target_arch="x86_64")]{match batch_ntt_simd_backend(a.len(),width){SimdBackend::Avx512=>{unsafe{ntt_simd::intt_batch_avx512::<_,false>(a,width)};return;}SimdBackend::Avx2=>{unsafe{ntt_simd::intt_batch_avx2::<_,false>(a,width)};return;}SimdBackend::Scalar=>{}}}intt_batch_scalar(a,width);}fn intt_batch_scalar(a:&mut[MInt],width:usize)where M:Montgomery32NttModulus{let n=a.len()/width;if n<=1{return;}let a=unsafe{std::slice::from_raw_parts_mut(a.as_mut_ptr().cast::(),a.len())};let mut v=1;let limit=if n.trailing_zeros()&1==1{n/2}else{n};let iimag=M::INFO.inv_root[2];while v(a0,a1);let a0na1=sub_scalar::(a0,a1);let a2pa3=add_scalar::(a2,a3);let a2na3iimag=mul_scalar::(sub_scalar::(a2,a3),iimag);*x0=add_scalar::(a0pa1,a2pa3);*x1=mul_scalar::(add_scalar::(a0na1,a2na3iimag),w1);*x2=mul_scalar::(sub_scalar::(a0pa1,a2pa3),w2);*x3=mul_scalar::(sub_scalar::(a0na1,a2na3iimag),w3);}let rate=&M::INFO.inv_rate3_packed[s.trailing_ones()as usize];w1=M::mod_mul(w1,rate[1]);w2=M::mod_mul(w2,rate[3]);w3=M::mod_mul(w3,rate[5]);}v<<=2;}if n.trailing_zeros()&1==1{let(l,r)=a.split_at_mut(n/2*width);for(x0,x1)in l.iter_mut().zip(r){let a0=*x0;let a1=*x1;*x0=add_scalar::(a0,a1);*x1=sub_scalar::(a0,a1);}}let inv=M::mod_inv(>::from(n as u32));for a in a{*a=M::mod_mul(*a,inv);}}fn ntt(a:&mut[MInt])where M:Montgomery32NttModulus{#[cfg(target_arch="x86_64")]match simd_backend(){SimdBackend::Avx512=>unsafe{ntt_simd::ntt_batch_avx512(a,1)},SimdBackend::Avx2=>unsafe{ntt_simd::ntt_batch_avx2::<_,true>(a,1)},SimdBackend::Scalar=>ntt_scalar(a),}#[cfg(not(target_arch="x86_64"))]ntt_scalar(a);}fn intt(a:&mut[MInt])where M:Montgomery32NttModulus{#[cfg(target_arch="x86_64")]match simd_backend(){SimdBackend::Avx512=>unsafe{ntt_simd::intt_batch_avx512::<_,true>(a,1)},SimdBackend::Avx2=>unsafe{ntt_simd::intt_batch_avx2::<_,true>(a,1)},SimdBackend::Scalar=>intt_scalar(a),}#[cfg(not(target_arch="x86_64"))]intt_scalar(a);}fn ntt_rows(a:&mut[MInt],width:usize)where M:Montgomery32NttModulus{for row in a.chunks_exact_mut(width){ntt(row);}}fn intt_rows(a:&mut[MInt],width:usize)where M:Montgomery32NttModulus{for row in a.chunks_exact_mut(width){intt(row);}}#[cfg(target_arch="x86_64")]fn use_block_ntt(len:usize)->bool where M:Montgomery32NttModulus{len>=64&&M::MOD(f:&mut[MInt],g:&[MInt])where M:Montgomery32NttModulus{assert!(f.len()<=g.len());crate::avx_helper!(@dispatch simd_backend,SimdBackend;unsafe{ntt_simd::pointwise_multiply_avx512(f,g)},unsafe{ntt_simd::pointwise_multiply_avx2(f,g)},{for(f,g)in f.iter_mut().zip(g.iter()){*f*=*g;}})}fn pointwise_multiply_add(sum:&mut[MInt],f:&[MInt],g:&[MInt])where M:Montgomery32NttModulus{crate::avx_helper!(@dispatch simd_backend,SimdBackend;unsafe{ntt_simd::pointwise_multiply_add_avx512(sum,f,g)},unsafe{ntt_simd::pointwise_multiply_add_avx2(sum,f,g)},{for((sum,f),g)in sum.iter_mut().zip(f.iter()).zip(g.iter()){*sum+=*f**g;}})}#[cfg(target_arch="x86_64")]#[allow(unsafe_op_in_unsafe_fn)]mod ntt_simd{#[cfg(target_arch="x86_64")]use super::*;#[cfg(target_arch="x86_64")]use std::arch::x86_64::*;#[inline]#[cfg(target_arch="x86_64")]fn normalize_scalar(x:u32)->u32 where M:Montgomery32NttModulus{if x>=M::MOD{x-M::MOD}else{x}}#[cfg(target_arch="x86_64")]mod convolution_avx2{#[cfg(target_arch="x86_64")]use super::*;#[inline]#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]unsafe fn load_block_avx2(a:*const u32,i:usize)->__m256i{_mm256_loadu_si256(a.add(i<<3)as*const __m256i)}#[inline]#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]unsafe fn store_block_avx2(a:*mut u32,i:usize,x:__m256i){_mm256_storeu_si256(a.add(i<<3)as*mut __m256i,x);}#[inline]#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]unsafe fn shrink_avx2(x:__m256i,modulus:__m256i)->__m256i{_mm256_min_epu32(x,_mm256_sub_epi32(x,modulus))}#[inline]#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]unsafe fn normalize_avx2(x:__m256i,modulus:__m256i,modulus2:__m256i)->__m256i{shrink_avx2(shrink_avx2(x,modulus2),modulus)}#[inline]#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]unsafe fn lazy_sub_avx2(x:__m256i,y:__m256i,modulus:__m256i)->__m256i{_mm256_add_epi32(x,_mm256_sub_epi32(modulus,y))}#[inline]#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]unsafe fn update_root_avx2(root:__m256i,rate:__m256i,modulus:__m256i)->__m256i{let correction=montgomery_simd::mul_epu32_256(root,rate);let product=_mm256_mul_epu32(root,_mm256_srli_epi64::<32>(rate));let correction=_mm256_mul_epu32(correction,modulus);shrink_avx2(_mm256_srli_epi64::<32>(_mm256_add_epi64(product,correction)),modulus)}#[inline]#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]unsafe fn packed_rate_avx2(index:usize,inverse:bool)->__m256i where M:Montgomery32NttModulus{let rate=if inverse{&M::INFO.inv_rate3_packed[index]}else{&M::INFO.rate3_packed[index]};_mm256_loadu_si256(rate.as_ptr().cast())}#[cfg(target_arch="x86_64")]macro_rules!block_butterfly{($forward:ident,$inverse:ident,$features:literal,$step:literal,$broadcast:expr,$load:ident,$store:ident,$add:ident,$sub:ident,$min:ident,$mul:path,$add_mod:path,$sub_mod:path)=>{#[inline]#[target_feature(enable=$features)]unsafe fn$forward(base:*mut u32,quarter:usize,root:__m256i)where M:Montgomery32NttModulus,{let modulus=($broadcast)(_mm256_set1_epi32(M::MOD as i32));let modulus2=($broadcast)(_mm256_set1_epi32((M::MOD*2)as i32));let imag=($broadcast)(_mm256_set1_epi32(M::INFO.root[2]as i32));let imag_r=($broadcast)(_mm256_set1_epi32(M::INFO.root[2].wrapping_mul(M::R)as i32));let indices=_mm256_setr_epi32(0,2,0,4,0,2,0,4);let r1=_mm256_permutevar8x32_epi32(root,indices);let r1_r=_mm256_permutevar8x32_epi32(montgomery_simd::mul_epu32_256(root,_mm256_set1_epi32(M::R as i32)),indices,);let r2=($broadcast)(_mm256_shuffle_epi32::<0x55>(r1));let nr3=($broadcast)(_mm256_shuffle_epi32::<0xff>(r1));let r2_r=($broadcast)(_mm256_shuffle_epi32::<0x55>(r1_r));let nr3_r=($broadcast)(_mm256_shuffle_epi32::<0xff>(r1_r));let r1=($broadcast)(r1);let r1_r=($broadcast)(r1_r);let shrink=|x,m|$min(x,$sub(x,m));let lazy_sub=|x,y,m|$add(x,$sub(m,y));let mut j=0;while j(base:*mut u32,quarter:usize,root:__m256i,inv:u32,group:usize,last:bool,)where M:Montgomery32NttModulus,{let modulus=($broadcast)(_mm256_set1_epi32(M::MOD as i32));let modulus2=($broadcast)(_mm256_set1_epi32((M::MOD*2)as i32));let imag=($broadcast)(_mm256_set1_epi32(M::INFO.root[2]as i32));let imag_r=($broadcast)(_mm256_set1_epi32(M::INFO.root[2].wrapping_mul(M::R)as i32));let indices=_mm256_setr_epi32(0,2,0,4,0,2,0,4);let r1=_mm256_permutevar8x32_epi32(root,indices);let r1_r=_mm256_permutevar8x32_epi32(montgomery_simd::mul_epu32_256(root,_mm256_set1_epi32(M::R as i32)),indices,);let r2=($broadcast)(_mm256_shuffle_epi32::<0x55>(r1));let r3=($broadcast)(_mm256_shuffle_epi32::<0xff>(r1));let r2_r=($broadcast)(_mm256_shuffle_epi32::<0x55>(r1_r));let r3_r=($broadcast)(_mm256_shuffle_epi32::<0xff>(r1_r));let r1=($broadcast)(r1);let r1_r=($broadcast)(r1_r);let shrink=|x,m|$min(x,$sub(x,m));let lazy_sub=|x,y,m|$add(x,$sub(m,y));let inv_vec=($broadcast)(_mm256_set1_epi32(inv as i32));let inv_r=($broadcast)(_mm256_set1_epi32(inv.wrapping_mul(M::R)as i32));let mut j=0;while j1{(shrink(h1,modulus2),$sub_mod(g0,g2,modulus2),$sub_mod(g1,g3,modulus2),)}else{($mul(h1,r1,r1_r,modulus),$mul(h2,r2,r2_r,modulus),$mul(h3,r3,r3_r,modulus),)};if last{y0=shrink(y0,modulus);y1=shrink(y1,modulus);y2=shrink(y2,modulus);y3=shrink(y3,modulus);}$store(base.add(p0).cast(),y0);$store(base.add(p0+(quarter<<3)).cast(),y1);$store(base.add(p0+(quarter<<4)).cast(),y2);$store(base.add(p0+quarter*24).cast(),y3);j+=$step;}}};}#[cfg(target_arch="x86_64")]block_butterfly!(ntt_butterfly_avx2,intt_butterfly_avx2,"avx2",1,|x|x,_mm256_loadu_si256,_mm256_storeu_si256,_mm256_add_epi32,_mm256_sub_epi32,_mm256_min_epu32,montgomery_simd::montgomery_mul_256_fixed,montgomery_simd::add_mod_256,montgomery_simd::sub_mod_256);#[cfg(target_arch="x86_64")]block_butterfly!(ntt_butterfly_avx512,intt_butterfly_avx512,"avx2,avx512f,avx512dq,avx512bw,avx512vl,avx512cd",2,_mm512_broadcast_i64x4,_mm512_loadu_si512,_mm512_storeu_si512,_mm512_add_epi32,_mm512_sub_epi32,_mm512_min_epu32,montgomery_simd::montgomery_mul_512_fixed,montgomery_simd::add_mod_512,montgomery_simd::sub_mod_512);#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]unsafe fn ntt_blocks_avx2(a:*mut u32,n:usize,input_len:usize)where M:Montgomery32NttModulus{let modulus=_mm256_set1_epi32(M::MOD as i32);let modulus2=_mm256_set1_epi32(M::MOD.wrapping_mul(2)as i32);let wide=matches!(simd_backend(),SimdBackend::Avx512);let imag=M::INFO.root[2];let imag_r=_mm256_set1_epi32(imag.wrapping_mul(M::R)as i32);let imag=_mm256_set1_epi32(imag as i32);let initial_root=_mm256_setr_epi32(M::INFO.root[3]as i32,0,M::INFO.root[2]as i32,0,(M::MOD-M::mod_mul(M::INFO.root[2],M::INFO.root[3]))as i32,0,0,0);let log_n=n.trailing_zeros()as usize;let mut roots=[initial_root;16];let nn=n>>(log_n&1);let tile_len=n.min(64);let mut stage_log=log_n&!1;let sparse=nn!=n&&input_len<=nn<<3;if nn!=n&&!sparse{let mut i=0;while i>2;let w1=M::INFO.root[3];let w3=M::MOD-M::mod_mul(w1,M::INFO.root[2]);let r1=_mm256_set1_epi32(w1 as i32);let r1_r=_mm256_set1_epi32(w1.wrapping_mul(M::R)as i32);let nr3=_mm256_set1_epi32(w3 as i32);let nr3_r=_mm256_set1_epi32(w3.wrapping_mul(M::R)as i32);let store=|offset,values:[__m256i;4]|{for(k,mut y)in values.into_iter().enumerate(){if quarter==1{y=normalize_avx2(y,modulus,modulus2);}store_block_avx2(a,offset+quarter*k,y);}};for j in 0..quarter{let x0=load_block_avx2(a,j);let x1=load_block_avx2(a,j+quarter);let x2=load_block_avx2(a,j+quarter*2);let x3=load_block_avx2(a,j+quarter*3);let h0=_mm256_add_epi32(x0,x2);let h1=_mm256_add_epi32(x1,x3);let h2=lazy_sub_avx2(x0,x2,modulus);let h3=montgomery_simd::montgomery_mul_256_fixed(lazy_sub_avx2(x1,x3,modulus),imag,imag_r,modulus);store(j,[montgomery_simd::add_mod_256(h0,h1,modulus2),lazy_sub_avx2(h0,h1,modulus2),_mm256_add_epi32(h2,h3),lazy_sub_avx2(h2,h3,modulus2)]);let g1=montgomery_simd::montgomery_mul_256_fixed(x1,r1,r1_r,modulus);let ng3=montgomery_simd::montgomery_mul_256_fixed(x3,nr3,nr3_r,modulus);let g2=montgomery_simd::montgomery_mul_256_fixed(x2,imag,imag_r,modulus);let h3=montgomery_simd::montgomery_mul_256_fixed(_mm256_add_epi32(g1,ng3),imag,imag_r,modulus);let h1=montgomery_simd::sub_mod_256(g1,ng3,modulus2);let h0=montgomery_simd::add_mod_256(x0,g2,modulus2);let h2=montgomery_simd::sub_mod_256(x0,g2,modulus2);store(nn+j,[_mm256_add_epi32(h0,h1),lazy_sub_avx2(h0,h1,modulus2),_mm256_add_epi32(h2,h3),lazy_sub_avx2(h2,h3,modulus2)]);}stage_log-=2;}if nn==n&&input_len<=n<<2{let quarter=n>>2;let mut i=0;while i=2{let quarter=1<<(stage_log-2);if wide&&quarter>1{ntt_butterfly_avx512::(a,quarter,initial_root);}else{ntt_butterfly_avx2::(a,quarter,initial_root);}stage_log-=2;}stage_log=max_stage_log.min(6);let mut tile=0;while tile>2;while quarter>1{let root_slot=(stage_log-2)>>1;let mut root=roots[root_slot];let mut i=if tile==0{group_len}else{0};let mut group=(tile+i)>>stage_log;while i((!group).trailing_zeros()as usize,false),modulus);let base=base.add(i<<3);if wide{ntt_butterfly_avx512::(base,quarter,current);}else{ntt_butterfly_avx2::(base,quarter,current);}i+=group_len;group+=1;}roots[root_slot]=root;group_len=quarter;quarter>>=2;stage_log-=2;}let r=_mm256_set1_epi32(M::R as i32);let store=|i,y|store_block_avx2(a,i,normalize_avx2(y,modulus,modulus2));let mut root=roots[0];let mut i=tile+if tile==0{4}else{0};while i((!(i>>2)).trailing_zeros()as usize,false),modulus);let r2=_mm256_shuffle_epi32::<0x55>(r1);let nr3=_mm256_shuffle_epi32::<0xff>(r1);let x0=load_block_avx2(a,i);let x1=load_block_avx2(a,i+1);let x2=load_block_avx2(a,i+2);let x3=load_block_avx2(a,i+3);let g1=montgomery_simd::montgomery_mul_256(x1,_mm256_shuffle_epi32::<0xa0>(r1),r,modulus);let ng3=montgomery_simd::montgomery_mul_256(x3,nr3,r,modulus);let g2=montgomery_simd::montgomery_mul_256(x2,r2,r,modulus);let g0=shrink_avx2(x0,modulus2);let h3=montgomery_simd::montgomery_mul_256_fixed(_mm256_add_epi32(g1,ng3),imag,imag_r,modulus);let h1=montgomery_simd::sub_mod_256(g1,ng3,modulus2);let h0=montgomery_simd::add_mod_256(g0,g2,modulus2);let h2=montgomery_simd::sub_mod_256(g0,g2,modulus2);store(i,_mm256_add_epi32(h0,h1));store(i+1,lazy_sub_avx2(h0,h1,modulus2));store(i+2,_mm256_add_epi32(h2,h3));store(i+3,lazy_sub_avx2(h2,h3,modulus2));i+=4;}roots[0]=root;tile+=tile_len;if tile(a:*mut u32,n:usize)where M:Montgomery32NttModulus{let inv=M::mod_inv(>::from(n as u32));let modulus=_mm256_set1_epi32(M::MOD as i32);let modulus2=_mm256_set1_epi32(M::MOD.wrapping_mul(2)as i32);let wide=matches!(simd_backend(),SimdBackend::Avx512);let initial_root=_mm256_setr_epi32(M::N1 as i32,0,M::N1 as i32,0,M::N1 as i32,0,0,0);let log_n=n.trailing_zeros()as usize;let mut roots=[initial_root;16];let nn=n>>(log_n&1);let tile_len=n.min(64);let inv_vec=_mm256_set1_epi32(inv as i32);roots[0]=inv_vec;let mut tile=0;while tile>1;let mut root=roots[root_slot];let mut group=(tile+i)>>stage_log;while i((!group).trailing_zeros()as usize,true),modulus);let base=base.add(i<<3);let last=group_len==n;if wide&&quarter>1{intt_butterfly_avx512::(base,quarter,current,inv,group,last);}else{intt_butterfly_avx2::(base,quarter,current,inv,group,last);}i+=group_len;group+=1;}roots[root_slot]=root;quarter=group_len;group_len<<=2;stage_log+=2;}tile+=tile_len;}if nn!=n{let mut i=0;while i__m256i{let even_m=_mm256_mul_epu32(even,r_vec);let odd_m=_mm256_mul_epu32(odd,r_vec);let even=_mm256_add_epi64(even,_mm256_mul_epu32(even_m,mod_vec));let odd=_mm256_add_epi64(odd,_mm256_mul_epu32(odd_m,mod_vec));_mm256_or_si256(_mm256_bsrli_epi128::<4>(even),odd)}#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]unsafe fn convolve_8_avx2(f:*mut u32,g:*const u32,n:usize)where M:Montgomery32NttModulus{let mod_vec=_mm256_set1_epi32(M::MOD as i32);let mod2_vec=_mm256_set1_epi32(M::MOD.wrapping_add(M::MOD)as i32);let r_vec=_mm256_set1_epi32(M::R as i32);let mut rr=M::N1;let mut i=0;while i::uninit();let work=work.as_mut_ptr().cast::();for(j,ww)in[rr,M::MOD.wrapping_mul(2)-rr,rr_i,M::MOD.wrapping_mul(2)-rr_i].into_iter().enumerate(){let k=i+j;let ff=load_block_avx2(f,k);let fw=montgomery_simd::montgomery_mul_256_fixed(ff,_mm256_set1_epi32(ww as i32),_mm256_set1_epi32(ww.wrapping_mul(M::R)as i32),mod_vec);_mm256_store_si256(work.add(j<<4).cast(),fw);_mm256_store_si256(work.add((j<<4)+8).cast(),ff);}let mut even=[_mm256_setzero_si256();4];let mut odd=[_mm256_setzero_si256();4];let mut l=0;while l<8{let mut j=0;while j<4{let x=_mm256_loadu_si256(work.add((j<<4)+8-l).cast());let y=_mm256_set1_epi32(*g.add(((i+j)<<3)+l)as i32);even[j]=_mm256_add_epi64(even[j],_mm256_mul_epu32(x,y));odd[j]=_mm256_add_epi64(odd[j],_mm256_mul_epu32(_mm256_bsrli_epi128::<4>(x),y));j+=1;}l+=1;}let mut j=0;while j<4{let x=reduce_sum_avx2(even[j],odd[j],r_vec,mod_vec);store_block_avx2(f,i+j,normalize_avx2(x,mod_vec,mod2_vec));j+=1;}i+=4;let rate=M::INFO.rate3_packed[(i>>2).trailing_zeros()as usize][1];rr=M::mod_mul(rr,rate);}}#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]pub unsafe fn transform_blocks_avx2(f:&mut[MInt],input_len:usize)where M:Montgomery32NttModulus{let n=f.len()>>3;let f=f.as_mut_ptr()as*mut u32;ntt_blocks_avx2::(f,n,input_len);}#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]pub unsafe fn multiply_blocks_avx2(f:&mut[MInt],g:Option<&[MInt]>)where M:Montgomery32NttModulus{let n=f.len()>>3;let f=f.as_mut_ptr()as*mut u32;let g=g.map_or(f.cast_const(),|g|g.as_ptr().cast());convolve_8_avx2::(f,g,n);}#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]pub unsafe fn inverse_transform_blocks_avx2(f:&mut[MInt])where M:Montgomery32NttModulus{let n=f.len()>>3;let f=f.as_mut_ptr()as*mut u32;intt_blocks_avx2::(f,n);}#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]pub unsafe fn convolve_blocks_avx2(f:&mut Vec>,g:&mut Vec>,same:bool,size:usize,input_len:(usize,usize))where M:Montgomery32NttModulus{let n=size>>3;let f_ptr=f.as_mut_ptr().cast();let g_ptr=if same{f_ptr}else{g.as_mut_ptr().cast()};ntt_blocks_avx2::(f_ptr,n,input_len.0);f.set_len(size);if!same{ntt_blocks_avx2::(g_ptr,n,input_len.1);}convolve_8_avx2::(f_ptr,g_ptr,n);intt_blocks_avx2::(f_ptr,n);}}#[cfg(target_arch="x86_64")]mod ntt_avx2{#[cfg(target_arch="x86_64")]use super::*;#[inline]#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]unsafe fn load_ntt_avx2(p:*const u32,step:usize)->__m256i{if step==4{_mm256_castsi128_si256(_mm_loadu_si128(p.cast()))}else{_mm256_loadu_si256(p.cast())}}#[inline]#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]unsafe fn store_ntt_avx2(p:*mut u32,value:__m256i,step:usize){if step==4{_mm_storeu_si128(p.cast(),_mm256_castsi256_si128(value));}else{_mm256_storeu_si256(p.cast(),value);}}#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]pub unsafe fn ntt_four_avx2(a:&mut[u32])where M:Montgomery32NttModulus{let roots=if INVERSE{&M::INFO.inv_root}else{&M::INFO.root};let rates=&const{let info=M::INFO;let(roots,rates)=if INVERSE{(info.inv_root,info.inv_rate3_packed)}else{(info.root,info.rate3_packed)};let mut result=[[0u32;4];32];let mut i=0;while i<32{let r=mod_mul(rates[i][1],roots[3],M::MOD,M::R);let r2=mod_mul(r,r,M::MOD,M::R);result[i]=[M::N1,r,r2,mod_mul(r2,r,M::MOD,M::R)];i+=1;}result};let one=_mm256_set1_epi32(M::N1 as i32);let modulus=_mm256_set1_epi32(M::MOD as i32);let modulus2=_mm256_set1_epi32((M::MOD*2)as i32);let r=_mm256_set1_epi32(M::R as i32);let root3=M::mod_mul(roots[2],roots[3])as i32;let step=_mm256_setr_epi32(M::N1 as i32,roots[3]as i32,roots[2]as i32,root3,M::N1 as i32,roots[3]as i32,roots[2]as i32,root3);let mut twiddle=_mm256_blend_epi32::<0xf0>(one,step);let imag=if INVERSE{_mm256_setr_epi32(M::N1 as i32,roots[2]as i32,M::N1 as i32,roots[2]as i32,M::N1 as i32,roots[2]as i32,M::N1 as i32,roots[2]as i32)}else{_mm256_setr_epi32(M::N1 as i32,M::N1 as i32,roots[2]as i32,roots[2]as i32,M::N1 as i32,M::N1 as i32,roots[2]as i32,roots[2]as i32)};for(s,a)in a.as_chunks_mut::<8>().0 .iter_mut().enumerate(){let mut x=_mm256_loadu_si256(a.as_ptr().cast());if!INVERSE{x=montgomery_simd::montgomery_mul_256(x,twiddle,r,modulus);}let pair=if INVERSE{let y=_mm256_shuffle_epi32::<0xb1>(x);let sum=montgomery_simd::add_mod_256(x,y,modulus2);let diff=montgomery_simd::sub_mod_256(x,y,modulus2);let sum=_mm256_shuffle_epi32::<0x88>(sum);let diff=_mm256_shuffle_epi32::<0x88>(diff);let diff=montgomery_simd::montgomery_mul_256(diff,imag,r,modulus);_mm256_unpacklo_epi64(sum,diff)}else{let y=_mm256_shuffle_epi32::<0x4e>(x);let sum=montgomery_simd::add_mod_256(x,y,modulus2);let diff=montgomery_simd::sub_mod_256(x,y,modulus2);_mm256_unpacklo_epi64(sum,diff)};let left=_mm256_shuffle_epi32::<0xa0>(pair);let mut right=_mm256_shuffle_epi32::<0xf5>(pair);if!INVERSE{right=montgomery_simd::montgomery_mul_256(right,imag,r,modulus);}let sum=montgomery_simd::add_mod_256(left,right,modulus2);let diff=montgomery_simd::sub_mod_256(left,right,modulus2);let mut value=_mm256_blend_epi32::<0xaa>(sum,diff);if INVERSE{value=_mm256_shuffle_epi32::<0xd8>(value);value=montgomery_simd::montgomery_mul_256(value,twiddle,r,modulus);}_mm256_storeu_si256(a.as_mut_ptr().cast(),value);let rate=_mm256_broadcastsi128_si256(_mm_loadu_si128(rates[s.trailing_ones()as usize+1].as_ptr().cast()));twiddle=montgomery_simd::montgomery_mul_256(twiddle,rate,r,modulus);}}#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]unsafe fn normalize_avx2(a:&mut[u32])where M:Montgomery32NttModulus{let mod_vec=_mm256_set1_epi32(M::MOD as i32);let mut i=0;while i+8<=a.len(){let x=_mm256_loadu_si256(a.as_ptr().add(i)as*const __m256i);let y=_mm256_min_epu32(x,_mm256_sub_epi32(x,mod_vec));_mm256_storeu_si256(a.as_mut_ptr().add(i)as*mut __m256i,y);i+=8;}while i(a[i]);i+=1;}}#[cfg(target_arch="x86_64")]pub unsafe fn add_vec_avx2(a:__m256i,b:__m256i,mod_vec:__m256i,mod2_vec:__m256i)->__m256i where M:Montgomery32NttModulus{if M::MOD(a:__m256i,b:__m256i,mod_vec:__m256i,mod2_vec:__m256i)->__m256i where M:Montgomery32NttModulus{if M::MOD(a:__m256i,b:__m256i,r_vec:__m256i,mod_vec:__m256i)->__m256i where M:Montgomery32NttModulus{if M::MOD(f:&mut[MInt],g:&[MInt])where M:Montgomery32NttModulus{let r_vec=_mm256_set1_epi32(M::R as i32);let mod_vec=_mm256_set1_epi32(M::MOD as i32);let mut i=0;while i+8<=f.len(){let a=_mm256_loadu_si256(f.as_ptr().add(i)as*const __m256i);let b=_mm256_loadu_si256(g.as_ptr().add(i)as*const __m256i);let x=montgomery_simd::montgomery_mul_256_canon(a,b,r_vec,mod_vec);_mm256_storeu_si256(f.as_mut_ptr().add(i)as*mut __m256i,x);i+=8;}while i(sum:&mut[MInt],f:&[MInt],g:&[MInt])where M:Montgomery32NttModulus{let r_vec=_mm256_set1_epi32(M::R as i32);let mod_vec=_mm256_set1_epi32(M::MOD as i32);let mut i=0;while i+8<=sum.len(){let s=_mm256_loadu_si256(sum.as_ptr().add(i).cast());let f=_mm256_loadu_si256(f.as_ptr().add(i).cast());let g=_mm256_loadu_si256(g.as_ptr().add(i).cast());let product=montgomery_simd::montgomery_mul_256_canon(f,g,r_vec,mod_vec);_mm256_storeu_si256(sum.as_mut_ptr().add(i).cast(),montgomery_simd::add_mod_256(s,product,mod_vec));i+=8;}while i(a:&mut[MInt],width:usize)where M:Montgomery32NttModulus{let n=a.len()/width;if n<=1{return;}let ptr=a.as_mut_ptr()as*mut u32;let a=std::slice::from_raw_parts_mut(ptr,a.len());let mod_vec=_mm256_set1_epi32(M::MOD as i32);let mod2_vec=_mm256_set1_epi32(M::MOD.wrapping_add(M::MOD)as i32);let r_vec=_mm256_set1_epi32(M::R as i32);let imag=M::INFO.root[2];let imag_vec=_mm256_set1_epi32(imag as i32);let mut v=n/2;if n.trailing_zeros()&1==1{let half=v*width;let step=if PARTIAL&&half==4{4}else{8};let mut i=0;while i+step<=half{let x0=load_ntt_avx2(a.as_ptr().add(i),step);let x1=load_ntt_avx2(a.as_ptr().add(half+i),step);let y0=add_vec_avx2::(x0,x1,mod_vec,mod2_vec);let y1=sub_vec_avx2::(x0,x1,mod_vec,mod2_vec);store_ntt_avx2(a.as_mut_ptr().add(i),y0,step);store_ntt_avx2(a.as_mut_ptr().add(half+i),y1,step);i+=step;}while i(x0,x1);a[half+i]=sub_scalar::(x0,x1);i+=1;}v>>=1;}while v>1{if width==1&&v==2&&a.len()>=8&&M::MOD(a);break;}let half=(v>>1)*width;let step=if PARTIAL&&half==4{4}else{8};let mut w1=M::N1;let mut w2=w1;let mut w3=w1;for(s,block)in a.chunks_exact_mut((v<<1)*width).enumerate(){let base=block.as_mut_ptr();let ll=base;let lr=base.add(half);let rl=base.add(v*width);let rr=base.add(v*width+half);let w1v=_mm256_set1_epi32(w1 as i32);let w2v=_mm256_set1_epi32(w2 as i32);let w3v=_mm256_set1_epi32(w3 as i32);let mut i=0;while i+step<=half{let x0=load_ntt_avx2(ll.add(i),step);let x1=load_ntt_avx2(lr.add(i),step);let x2=load_ntt_avx2(rl.add(i),step);let x3=load_ntt_avx2(rr.add(i),step);let(a1,a2,a3)=if s==0{(x1,x2,x3)}else{(mul_vec_avx2::(x1,w1v,r_vec,mod_vec),mul_vec_avx2::(x2,w2v,r_vec,mod_vec),mul_vec_avx2::(x3,w3v,r_vec,mod_vec))};let a0pa2=add_vec_avx2::(x0,a2,mod_vec,mod2_vec);let a0na2=sub_vec_avx2::(x0,a2,mod_vec,mod2_vec);let a1pa3=add_vec_avx2::(a1,a3,mod_vec,mod2_vec);let a1na3=sub_vec_avx2::(a1,a3,mod_vec,mod2_vec);let a1na3imag=mul_vec_avx2::(a1na3,imag_vec,r_vec,mod_vec);let y0=add_vec_avx2::(a0pa2,a1pa3,mod_vec,mod2_vec);let y1=sub_vec_avx2::(a0pa2,a1pa3,mod_vec,mod2_vec);let y2=add_vec_avx2::(a0na2,a1na3imag,mod_vec,mod2_vec);let y3=sub_vec_avx2::(a0na2,a1na3imag,mod_vec,mod2_vec);store_ntt_avx2(ll.add(i),y0,step);store_ntt_avx2(lr.add(i),y1,step);store_ntt_avx2(rl.add(i),y2,step);store_ntt_avx2(rr.add(i),y3,step);i+=step;}while i(*lr.add(i),w1);let a2=mul_scalar::(*rl.add(i),w2);let a3=mul_scalar::(*rr.add(i),w3);let a0pa2=add_scalar::(a0,a2);let a0na2=sub_scalar::(a0,a2);let a1pa3=add_scalar::(a1,a3);let a1na3=sub_scalar::(a1,a3);let a1na3imag=mul_scalar::(a1na3,imag);*ll.add(i)=add_scalar::(a0pa2,a1pa3);*lr.add(i)=sub_scalar::(a0pa2,a1pa3);*rl.add(i)=add_scalar::(a0na2,a1na3imag);*rr.add(i)=sub_scalar::(a0na2,a1na3imag);i+=1;}let rate=&M::INFO.rate3_packed[s.trailing_ones()as usize];w1=M::mod_mul(w1,rate[1]);w2=M::mod_mul(w2,rate[3]);w3=M::mod_mul(w3,rate[5]);}v>>=2;}normalize_avx2::(a);}#[inline]#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]pub unsafe fn intt_batch_avx2(a:&mut[MInt],width:usize)where M:Montgomery32NttModulus{let n=a.len()/width;if n<=1{return;}let ptr=a.as_mut_ptr()as*mut u32;let a=std::slice::from_raw_parts_mut(ptr,a.len());let mod_vec=_mm256_set1_epi32(M::MOD as i32);let mod2_vec=_mm256_set1_epi32(M::MOD.wrapping_add(M::MOD)as i32);let r_vec=_mm256_set1_epi32(M::R as i32);let iimag=M::INFO.inv_root[2];let iimag_vec=_mm256_set1_epi32(iimag as i32);let mut v=1;if width==1&&a.len()>=8&&M::MOD(a);v=4;}let limit=if n.trailing_zeros()&1==1{n/2}else{n};while v(x0,x1,mod_vec,mod2_vec);let a0na1=sub_vec_avx2::(x0,x1,mod_vec,mod2_vec);let a2pa3=add_vec_avx2::(x2,x3,mod_vec,mod2_vec);let a2na3=sub_vec_avx2::(x2,x3,mod_vec,mod2_vec);let a2na3iimag=mul_vec_avx2::(a2na3,iimag_vec,r_vec,mod_vec);let y0=add_vec_avx2::(a0pa1,a2pa3,mod_vec,mod2_vec);let y1=add_vec_avx2::(a0na1,a2na3iimag,mod_vec,mod2_vec);let y2=sub_vec_avx2::(a0pa1,a2pa3,mod_vec,mod2_vec);let y3=sub_vec_avx2::(a0na1,a2na3iimag,mod_vec,mod2_vec);let(y1,y2,y3)=if s==0{(y1,y2,y3)}else{(mul_vec_avx2::(y1,w1v,r_vec,mod_vec),mul_vec_avx2::(y2,w2v,r_vec,mod_vec),mul_vec_avx2::(y3,w3v,r_vec,mod_vec))};store_ntt_avx2(ll.add(i),y0,step);store_ntt_avx2(lr.add(i),y1,step);store_ntt_avx2(rl.add(i),y2,step);store_ntt_avx2(rr.add(i),y3,step);i+=step;}while i(a0,a1);let a0na1=sub_scalar::(a0,a1);let a2pa3=add_scalar::(a2,a3);let a2na3iimag=mul_scalar::(sub_scalar::(a2,a3),iimag);*ll.add(i)=add_scalar::(a0pa1,a2pa3);*lr.add(i)=mul_scalar::(add_scalar::(a0na1,a2na3iimag),w1);*rl.add(i)=mul_scalar::(sub_scalar::(a0pa1,a2pa3),w2);*rr.add(i)=mul_scalar::(sub_scalar::(a0na1,a2na3iimag),w3);i+=1;}let rate=&M::INFO.inv_rate3_packed[s.trailing_ones()as usize];w1=M::mod_mul(w1,rate[1]);w2=M::mod_mul(w2,rate[3]);w3=M::mod_mul(w3,rate[5]);}v<<=2;}if n.trailing_zeros()&1==1{let half=(n>>1)*width;let step=if PARTIAL&&half==4{4}else{8};let mut i=0;while i+step<=half{let x0=load_ntt_avx2(a.as_ptr().add(i),step);let x1=load_ntt_avx2(a.as_ptr().add(half+i),step);let y0=add_vec_avx2::(x0,x1,mod_vec,mod2_vec);let y1=sub_vec_avx2::(x0,x1,mod_vec,mod2_vec);store_ntt_avx2(a.as_mut_ptr().add(i),y0,step);store_ntt_avx2(a.as_mut_ptr().add(half+i),y1,step);i+=step;}while i(x0,x1);a[half+i]=sub_scalar::(x0,x1);i+=1;}}let inv=M::mod_inv(>::from(n as u32));let inv_vec=_mm256_set1_epi32(inv as i32);let mut i=0;while i+8<=a.len(){let x=_mm256_loadu_si256(a.as_ptr().add(i)as*const __m256i);let y=montgomery_simd::montgomery_mul_256_canon(x,inv_vec,r_vec,mod_vec);_mm256_storeu_si256(a.as_mut_ptr().add(i)as*mut __m256i,y);i+=8;}while i(a:&mut[u32])where M:Montgomery32NttModulus{let mod_vec=_mm512_set1_epi32(M::MOD as i32);let mut i=0;while i+16<=a.len(){let x=_mm512_loadu_si512(a.as_ptr().add(i)as*const __m512i);let y=_mm512_min_epu32(x,_mm512_sub_epi32(x,mod_vec));_mm512_storeu_si512(a.as_mut_ptr().add(i)as*mut __m512i,y);i+=16;}while i(a[i]);i+=1;}}#[cfg(target_arch="x86_64")]unsafe fn add_vec_avx512(a:__m512i,b:__m512i,mod_vec:__m512i,mod2_vec:__m512i)->__m512i where M:Montgomery32NttModulus{if M::MOD(a:__m512i,b:__m512i,mod_vec:__m512i,mod2_vec:__m512i)->__m512i where M:Montgomery32NttModulus{if M::MOD(a:__m512i,b:__m512i,r_vec:__m512i,mod_vec:__m512i)->__m512i where M:Montgomery32NttModulus{if M::MOD(f:&mut[MInt],g:&[MInt])where M:Montgomery32NttModulus{let r_vec=_mm512_set1_epi32(M::R as i32);let mod_vec=_mm512_set1_epi32(M::MOD as i32);let mut i=0;while i+16<=f.len(){let a=_mm512_loadu_si512(f.as_ptr().add(i)as*const __m512i);let b=_mm512_loadu_si512(g.as_ptr().add(i)as*const __m512i);let x=montgomery_simd::montgomery_mul_512_canon(a,b,r_vec,mod_vec);_mm512_storeu_si512(f.as_mut_ptr().add(i)as*mut __m512i,x);i+=16;}while i(sum:&mut[MInt],f:&[MInt],g:&[MInt])where M:Montgomery32NttModulus{let r_vec=_mm512_set1_epi32(M::R as i32);let mod_vec=_mm512_set1_epi32(M::MOD as i32);let mut i=0;while i+16<=sum.len(){let s=_mm512_loadu_si512(sum.as_ptr().add(i).cast());let f=_mm512_loadu_si512(f.as_ptr().add(i).cast());let g=_mm512_loadu_si512(g.as_ptr().add(i).cast());let product=montgomery_simd::montgomery_mul_512_canon(f,g,r_vec,mod_vec);_mm512_storeu_si512(sum.as_mut_ptr().add(i).cast(),montgomery_simd::add_mod_512(s,product,mod_vec));i+=16;}while i(a:&mut[MInt],width:usize)where M:Montgomery32NttModulus{let n=a.len()/width;if n<=1{return;}let ptr=a.as_mut_ptr()as*mut u32;let a=std::slice::from_raw_parts_mut(ptr,a.len());let mod_vec=_mm512_set1_epi32(M::MOD as i32);let mod2_vec=_mm512_set1_epi32(M::MOD.wrapping_add(M::MOD)as i32);let r_vec=_mm512_set1_epi32(M::R as i32);let imag=M::INFO.root[2];let imag_vec=_mm512_set1_epi32(imag as i32);let mut v=n/2;if n.trailing_zeros()&1==1{let half=v*width;let mut i=0;while i+16<=half{let x0=_mm512_loadu_si512(a.as_ptr().add(i)as*const __m512i);let x1=_mm512_loadu_si512(a.as_ptr().add(half+i)as*const __m512i);let y0=add_vec_avx512::(x0,x1,mod_vec,mod2_vec);let y1=sub_vec_avx512::(x0,x1,mod_vec,mod2_vec);_mm512_storeu_si512(a.as_mut_ptr().add(i)as*mut __m512i,y0);_mm512_storeu_si512(a.as_mut_ptr().add(half+i)as*mut __m512i,y1);i+=16;}while i>=1;}while v>1{if width==1&&v==2&&a.len()>=8&&M::MOD(a);break;}let half=(v>>1)*width;let mut w1=M::N1;for(s,block)in a.chunks_exact_mut((v<<1)*width).enumerate(){let base=block.as_mut_ptr();let ll=base;let lr=base.add(half);let rl=base.add(v*width);let rr=base.add(v*width+half);let w2=M::mod_mul(w1,w1);let w3=M::mod_mul(w2,w1);let w1v=_mm512_set1_epi32(w1 as i32);let w2v=_mm512_set1_epi32(w2 as i32);let w3v=_mm512_set1_epi32(w3 as i32);let mut i=0;while i+16<=half{let x0=_mm512_loadu_si512(ll.add(i)as*const __m512i);let x1=_mm512_loadu_si512(lr.add(i)as*const __m512i);let x2=_mm512_loadu_si512(rl.add(i)as*const __m512i);let x3=_mm512_loadu_si512(rr.add(i)as*const __m512i);let(a1,a2,a3)=if s==0{(x1,x2,x3)}else{(mul_vec_avx512::(x1,w1v,r_vec,mod_vec),mul_vec_avx512::(x2,w2v,r_vec,mod_vec),mul_vec_avx512::(x3,w3v,r_vec,mod_vec))};let a0pa2=add_vec_avx512::(x0,a2,mod_vec,mod2_vec);let a0na2=sub_vec_avx512::(x0,a2,mod_vec,mod2_vec);let a1pa3=add_vec_avx512::(a1,a3,mod_vec,mod2_vec);let a1na3=sub_vec_avx512::(a1,a3,mod_vec,mod2_vec);let a1na3imag=mul_vec_avx512::(a1na3,imag_vec,r_vec,mod_vec);let y0=add_vec_avx512::(a0pa2,a1pa3,mod_vec,mod2_vec);let y1=sub_vec_avx512::(a0pa2,a1pa3,mod_vec,mod2_vec);let y2=add_vec_avx512::(a0na2,a1na3imag,mod_vec,mod2_vec);let y3=sub_vec_avx512::(a0na2,a1na3imag,mod_vec,mod2_vec);_mm512_storeu_si512(ll.add(i)as*mut __m512i,y0);_mm512_storeu_si512(lr.add(i)as*mut __m512i,y1);_mm512_storeu_si512(rl.add(i)as*mut __m512i,y2);_mm512_storeu_si512(rr.add(i)as*mut __m512i,y3);i+=16;}while i(*ll.add(i));let a1=M::mod_mul(normalize_scalar::(*lr.add(i)),w1);let a2=M::mod_mul(normalize_scalar::(*rl.add(i)),w2);let a3=M::mod_mul(normalize_scalar::(*rr.add(i)),w3);let a0pa2=M::mod_add(a0,a2);let a0na2=M::mod_sub(a0,a2);let a1pa3=M::mod_add(a1,a3);let a1na3=M::mod_sub(a1,a3);let a1na3imag=M::mod_mul(a1na3,imag);*ll.add(i)=M::mod_add(a0pa2,a1pa3);*lr.add(i)=M::mod_sub(a0pa2,a1pa3);*rl.add(i)=M::mod_add(a0na2,a1na3imag);*rr.add(i)=M::mod_sub(a0na2,a1na3imag);i+=1;}w1=M::mod_mul(w1,M::INFO.rate3_packed[s.trailing_ones()as usize][1]);}v>>=2;}normalize_avx512::(a);}#[inline]#[target_feature(enable="avx512f,avx512dq,avx512cd,avx512bw,avx512vl")]#[cfg(target_arch="x86_64")]pub unsafe fn intt_batch_avx512(a:&mut[MInt],width:usize)where M:Montgomery32NttModulus{let width=if SINGLE{1}else{width};let n=a.len()/width;if n<=1{return;}let ptr=a.as_mut_ptr()as*mut u32;let a=std::slice::from_raw_parts_mut(ptr,a.len());let mod_vec=_mm512_set1_epi32(M::MOD as i32);let mod2_vec=_mm512_set1_epi32(M::MOD.wrapping_add(M::MOD)as i32);let r_vec=_mm512_set1_epi32(M::R as i32);let iimag=M::INFO.inv_root[2];let iimag_vec=_mm512_set1_epi32(iimag as i32);let mut v=1;if width==1&&a.len()>=8&&M::MOD(a);v=4;}let limit=if n.trailing_zeros()&1==1{n/2}else{n};while v(x0,x1,mod_vec,mod2_vec);let a0na1=sub_vec_avx512::(x0,x1,mod_vec,mod2_vec);let a2pa3=add_vec_avx512::(x2,x3,mod_vec,mod2_vec);let a2na3=sub_vec_avx512::(x2,x3,mod_vec,mod2_vec);let a2na3iimag=mul_vec_avx512::(a2na3,iimag_vec,r_vec,mod_vec);let y0=add_vec_avx512::(a0pa1,a2pa3,mod_vec,mod2_vec);let y1=add_vec_avx512::(a0na1,a2na3iimag,mod_vec,mod2_vec);let y2=sub_vec_avx512::(a0pa1,a2pa3,mod_vec,mod2_vec);let y3=sub_vec_avx512::(a0na1,a2na3iimag,mod_vec,mod2_vec);let(y1,y2,y3)=if s==0{(y1,y2,y3)}else{(mul_vec_avx512::(y1,w1v,r_vec,mod_vec),mul_vec_avx512::(y2,w2v,r_vec,mod_vec),mul_vec_avx512::(y3,w3v,r_vec,mod_vec))};_mm512_storeu_si512(ll.add(i)as*mut __m512i,y0);_mm512_storeu_si512(lr.add(i)as*mut __m512i,y1);_mm512_storeu_si512(rl.add(i)as*mut __m512i,y2);_mm512_storeu_si512(rr.add(i)as*mut __m512i,y3);i+=16;}while i(a0,a1);let a0na1=sub_scalar::(a0,a1);let a2pa3=add_scalar::(a2,a3);let a2na3iimag=mul_scalar::(sub_scalar::(a2,a3),iimag);*ll.add(i)=add_scalar::(a0pa1,a2pa3);*lr.add(i)=mul_scalar::(add_scalar::(a0na1,a2na3iimag),w1);*rl.add(i)=mul_scalar::(sub_scalar::(a0pa1,a2pa3),w2);*rr.add(i)=mul_scalar::(sub_scalar::(a0na1,a2na3iimag),w3);i+=1;}let rate=&M::INFO.inv_rate3_packed[s.trailing_ones()as usize];w1=M::mod_mul(w1,rate[1]);w2=M::mod_mul(w2,rate[3]);w3=M::mod_mul(w3,rate[5]);}v<<=2;}if n.trailing_zeros()&1==1{let half=n/2*width;let mut i=0;while i+16<=half{let x0=_mm512_loadu_si512(a.as_ptr().add(i)as*const __m512i);let x1=_mm512_loadu_si512(a.as_ptr().add(half+i)as*const __m512i);let y0=add_vec_avx512::(x0,x1,mod_vec,mod2_vec);let y1=sub_vec_avx512::(x0,x1,mod_vec,mod2_vec);_mm512_storeu_si512(a.as_mut_ptr().add(i)as*mut __m512i,y0);_mm512_storeu_si512(a.as_mut_ptr().add(half+i)as*mut __m512i,y1);i+=16;}while i(x0,x1);a[half+i]=sub_scalar::(x0,x1);i+=1;}}let inv=M::mod_inv(>::from(n as u32));let inv_vec=_mm512_set1_epi32(inv as i32);let mut i=0;while i+16<=a.len(){let x=_mm512_loadu_si512(a.as_ptr().add(i)as*const __m512i);let y=montgomery_simd::montgomery_mul_512_canon(x,inv_vec,r_vec,mod_vec);_mm512_storeu_si512(a.as_mut_ptr().add(i)as*mut __m512i,y);i+=16;}while i(a:&[T],b:&[T])->Vecwhere T:Copy+Zero+AddAssign+Mul{if a.is_empty()&&b.is_empty(){return Vec::new();}let len=a.len()+b.len()-1;let mut c=vec![T::zero();len];let(a,b)=if a.len()>=b.len(){(a,b)}else{(b,a)};if!BLOCKED{for(i,&b)in b.iter().enumerate(){for(a,c)in a.iter().zip(&mut c[i..i+a.len()]).rev(){*c+=*a*b;}}return c;}let reverse=cfg!(target_arch="x86_64")&&(a.as_ptr().addr()^c.as_ptr().addr())&4095==0;for(block,a)in a.chunks(1024).enumerate(){for(i,&b)in b.iter().enumerate(){let start=block*1024+i;let pairs=a.iter().zip(&mut c[start..start+a.len()]);if reverse{for(a,c)in pairs.rev(){*c+=*a*b;}}else{for(a,c)in pairs{*c+=*a*b;}}}}c}fn convolve_karatsuba(a:&[T],b:&[T])->Vecwhere T:Copy+Zero+AddAssign+SubAssign+Mul{if a.len().min(b.len())<=30{return convolve_naive::<_,{!cfg!(target_arch="x86_64")}>(a,b);}let block_len=a.len().min(b.len()).next_power_of_two();if a.len().max(b.len())>block_len*4{let(a,b)=if a.len()>=b.len(){(a,b)}else{(b,a)};let mut result=vec![T::zero();a.len()+b.len()-1];for(i,a)in a.chunks(block_len).enumerate(){for(value,product)in result[i*block_len..].iter_mut().zip(convolve_karatsuba(a,b)){*value+=product;}}return result;}let m=a.len().max(b.len()).div_ceil(2);let(a0,a1)=if a.len()<=m{(a,&[][..])}else{a.split_at(m)};let(b0,b1)=if b.len()<=m{(b,&[][..])}else{b.split_at(m)};let f00=convolve_karatsuba(a0,b0);let f11=convolve_karatsuba(a1,b1);let mut a0a1=a0.to_vec();for(a0a1,&a1)in a0a1.iter_mut().zip(a1){*a0a1+=a1;}let mut b0b1=b0.to_vec();for(b0b1,&b1)in b0b1.iter_mut().zip(b1){*b0b1+=b1;}let mut f01=convolve_karatsuba(&a0a1,&b0b1);for(f01,&f00)in f01.iter_mut().zip(&f00){*f01-=f00;}for(f01,&f11)in f01.iter_mut().zip(&f11){*f01-=f11;}let mut c=vec![T::zero();a.len()+b.len()-1];for(c,&f00)in c.iter_mut().zip(&f00){*c+=f00;}for(c,&f01)in c[m..].iter_mut().zip(&f01){*c+=f01;}for(c,&f11)in c[m<<1..].iter_mut().zip(&f11){*c+=f11;}c}#[cold]fn convolve_large_ntt(a:Vec>,b:Vec>)->Vec>where M:Montgomery32NttModulus{let len=a.len()+b.len()-1;let ntt_len=1usize<]|{let mut f=Vec::with_capacity(ntt_len);advise_huge_pages(&mut f);f.extend_from_slice(a);Convolve::::transform_ntt(f,ntt_len)};let fa:Vec<_>=a.chunks(block_len).map(transform).collect();let fb:Option>=if same{None}else{Some(b.chunks(block_len).map(transform).collect())};let b_blocks=fb.as_ref().map_or(fa.len(),Vec::len);let mut result=vec![MInt::::zero();len];for diagonal in 0..fa.len()+b_blocks-1{let mut spectrum=vec![MInt::::zero();ntt_len];let start=diagonal.saturating_sub(b_blocks-1);for i in start..=diagonal.min(fa.len()-1){let j=diagonal-i;let g=if let Some(fb)=&fb{&fb[j]}else{&fa[j]};pointwise_multiply_add(&mut spectrum,&fa[i],g);}spectrum=Convolve::::inverse_transform_ntt(spectrum,ntt_len);let offset=diagonal*block_len;for(result,value)in result[offset..].iter_mut().zip(spectrum){*result+=value;}}result}implConvolveSteps for Convolvewhere M:Montgomery32NttModulus{const CYCLIC:bool=true;type T=Vec>;type F=Vec>;fn length(t:&Self::T)->usize{t.len()}fn transform(mut t:Self::T,len:usize)->Self::F{#[cfg(target_arch="x86_64")]let input_len=t.len();t.resize_with(len.max(1).next_power_of_two(),Zero::zero);#[cfg(target_arch="x86_64")]if use_block_ntt::(t.len()){unsafe{ntt_simd::transform_blocks_avx2(&mut t,input_len)};return t;}ntt(&mut t);t}fn inverse_transform(mut f:Self::F,len:usize)->Self::T{#[cfg(target_arch="x86_64")]if use_block_ntt::(f.len()){unsafe{ntt_simd::inverse_transform_blocks_avx2(&mut f)};f.truncate(len);return f;}intt(&mut f);f.truncate(len);f}fn multiply(f:&mut Self::F,g:&Self::F){assert_eq!(f.len(),g.len());#[cfg(target_arch="x86_64")]if use_block_ntt::(f.len()){unsafe{ntt_simd::multiply_blocks_avx2(f,Some(g))};return;}pointwise_multiply(f,g);}fn square(t:Self::T,len:usize)->Self::T{let mut f=Self::transform(t,len);#[cfg(target_arch="x86_64")]if use_block_ntt::(f.len()){unsafe{ntt_simd::multiply_blocks_avx2(&mut f,None)};return Self::inverse_transform(f,len);}let g=f.clone();Self::multiply(&mut f,&g);Self::inverse_transform(f,len)}fn convolve(mut a:Self::T,mut b:Self::T)->Self::T{let len=(Self::length(&a)+Self::length(&b)).saturating_sub(1);let(threshold,naive_threshold)=(100,60);#[cfg(target_arch="x86_64")]let(threshold,naive_threshold)=if use_block_ntt::(64){(60,if M::RANK>=13&&len<=8192{18}else if M::RANK>=19&&len<=262144{32}else{34})}else{(threshold,naive_threshold)};if Self::length(&a).max(Self::length(&b))<=threshold{return convolve_karatsuba(&a,&b);}if Self::length(&a).min(Self::length(&b))<=naive_threshold{return convolve_naive::<_,true>(&a,&b);}let size=len.max(1).next_power_of_two();let max_size=1usize<(size){max_size<<3}else{max_size};if size>max_size{return convolve_large_ntt(a,b);}if len<=size/2+2{let xa=a.pop().unwrap();let xb=b.pop().unwrap();let mut c=vec![MInt::::zero();len];*c.last_mut().unwrap()=xa*xb;for(a,c)in a.iter().zip(&mut c[b.len()..]){*c+=*a*xb;}for(b,c)in b.iter().zip(&mut c[a.len()..]){*c+=*b*xa;}let d=Self::convolve(a,b);for(d,c)in d.into_iter().zip(&mut c){*c+=d;}return c;}let same=a==b;#[cfg(target_arch="x86_64")]if use_block_ntt::(size){let input_len=(a.len(),b.len());for x in[&mut a,&mut b]{x.reserve(size-x.len());advise_huge_pages(x);let initialized=if x.len()<=size/2{size/2}else{size};x.resize_with(initialized,Zero::zero);}unsafe{ntt_simd::convolve_blocks_avx2(&mut a,&mut b,same,size,input_len)};a.truncate(len);return a;}let mut a=Self::transform(a,len);if same{for a in a.iter_mut(){*a*=*a;}}else{let b=Self::transform(b,len);Self::multiply(&mut a,&b);}Self::inverse_transform(a,len)}}type MVec=Vec>;fn convert_crt_input(t:MVec,capacity:usize)->(MVec,MVec,MVec)where M:MIntConvert,N1:Montgomery32NttModulus,N2:Montgomery32NttModulus,N3:Montgomery32NttModulus{let mut f=(MVec::::with_capacity(capacity),MVec::::with_capacity(capacity),MVec::::with_capacity(capacity));advise_huge_pages(&mut f.0);advise_huge_pages(&mut f.1);advise_huge_pages(&mut f.2);for t in t{let t:u32=t.into();f.0 .push(t.into());f.1 .push(t.into());f.2 .push(t.into());}f}fn reconstruct_mint_crt(f:(MVec,MVec,MVec))->MVecwhere M:MIntConvert+MIntConvert,N1:Montgomery32NttModulus,N2:Montgomery32NttModulus,N3:Montgomery32NttModulus{let t1=MInt::::new(N1::get_mod()).inv();let m1_3=MInt::::new(N1::get_mod());let t2=(m1_3*MInt::::new(N2::get_mod())).inv();let modulus=>::mod_into()as u64;let m1=N1::get_mod()as u64;let m2=m1*N2::get_mod()as u64%modulus;let fits_u64=(N1::get_mod()-1)as u128+(N2::get_mod()-1)as u128*m1 as u128+(N3::get_mod()-1)as u128*m2 as u128<=u64::MAX as u128;f.0 .into_iter().zip(f.1).zip(f.2).map(|((c1,c2),c3)|{let d1=c1.inner();let d2=((c2-MInt::::from(d1))*t1).inner();let x=MInt::::new(d1)+MInt::::new(d2)*m1_3;let d3=((c3-x)*t2).inner();let value=if fits_u64{(d1 as u64+d2 as u64*m1+d3 as u64*m2)%modulus}else{((d1 as u128+d2 as u128*m1 as u128+d3 as u128*m2 as u128)%modulus as u128)as u64};MInt::::from(value as u32)}).collect()}implConvolveSteps for Convolve<(M,(N1,N2,N3))>where M:MIntConvert+MIntConvert,N1:Montgomery32NttModulus,N2:Montgomery32NttModulus,N3:Montgomery32NttModulus{type T=MVec;type F=(MVec,MVec,MVec);fn length(t:&Self::T)->usize{t.len()}fn transform(t:Self::T,len:usize)->Self::F{let npot=len.max(1).next_power_of_two();let f=convert_crt_input(t,npot);(Convolve::::transform(f.0,npot),Convolve::::transform(f.1,npot),Convolve::::transform(f.2,npot))}fn inverse_transform(f:Self::F,len:usize)->Self::T{reconstruct_mint_crt((Convolve::::inverse_transform(f.0,len),Convolve::::inverse_transform(f.1,len),Convolve::::inverse_transform(f.2,len)))}fn multiply(f:&mut Self::F,g:&Self::F){Convolve::::multiply(&mut f.0,&g.0);Convolve::::multiply(&mut f.1,&g.1);Convolve::::multiply(&mut f.2,&g.2);}fn convolve(a:Self::T,b:Self::T)->Self::T{let max_len=Self::length(&a).max(Self::length(&b));let min_len=Self::length(&a).min(Self::length(&b));let(balanced,short)=crate::avx_helper!(@dispatch_avx2_fma(30,10),(384,128));if max_len<=balanced||min_len<=short{return convolve_karatsuba(&a,&b);}let fft_limit=crate::avx_helper!(@dispatch_avx2_fma 1usize<<((1u64<<50)/>::mod_into()as u64).ilog2().min(20),0);let convolve=|a:Self::T,b:Self::T|{let fft_len=(a.len()+b.len()-1).next_power_of_two();if fft_len<=256&&a.len()*b.len()<=fft_len*8{return convolve_karatsuba(&a,&b);}if fft_len<=fft_limit{crate::avx_helper!(@dispatch_avx2_fma return unsafe{convolve_mint_avx2(a,b)},());}convolve_mint_crt::(a,b)};let block_len=min_len.next_power_of_two()*8-min_len+1;let block_len=if min_len<=fft_limit/2{block_len.min(fft_limit-min_len+1)}else{block_len};if max_len<=block_len{return convolve(a,b);}let(a,b)=if a.len()>=b.len(){(a,b)}else{(b,a)};let mut result=vec![MInt::::zero();a.len()+b.len()-1];for(i,a)in a.chunks(block_len).enumerate(){let product=convolve(a.to_vec(),b.clone());for(value,product)in result[i*block_len..].iter_mut().zip(product){*value+=product;}}result}}fn convolve_mint_crt(a:MVec,b:MVec)->MVecwhere M:MIntConvert+MIntConvert,N1:Montgomery32NttModulus,N2:Montgomery32NttModulus,N3:Montgomery32NttModulus{let convolve=|a:MVec,b:MVec|{let a_len=a.len();let b_len=b.len();let a=convert_crt_input(a,a_len);let b=convert_crt_input(b,b_len);reconstruct_mint_crt((Convolve::::convolve(a.0,b.0),Convolve::::convolve(a.1,b.1),Convolve::::convolve(a.2,b.2)))};let modulus=>::mod_into()as u128;let capacity=N1::MOD as u128*N2::MOD as u128*N3::MOD as u128;if a.len().min(b.len())as u128*(modulus-1).pow(2)(&a,&b);}let mut result=vec![MInt::::zero();a.len()+b.len()-1];for(i,a)in a.chunks(block_len).enumerate(){for(j,b)in b.chunks(block_len).enumerate(){let product=convolve(a.to_vec(),b.to_vec());for(value,product)in result[(i+j)*block_len..].iter_mut().zip(product){*value+=product;}}}result}implConvolveSteps for Convolve<(u64,(N1,N2,N3))>where N1:Montgomery32NttModulus,N2:Montgomery32NttModulus,N3:Montgomery32NttModulus{type T=Vec;type F=([MVec;3],[MVec;3],[MVec;3]);fn length(t:&Self::T)->usize{t.len()}fn transform(t:Self::T,len:usize)->Self::F{let npot=len.max(1).next_power_of_two();assert!(npot<=1usize<(t:&[u64],len:usize,bits:usize,parts:usize)->[MVec;3]{std::array::from_fn(|part|{if part>=parts{return Vec::new();}Convolve::::transform(t.iter().map(|&t|MInt::from((t>>(part*bits))&((1u64<Self::T{let bits=if f.0[2].is_empty(){32}else{22};let t1=MInt::::new(N1::get_mod()).inv();let m1=N1::get_mod()as u64;let m1_3=MInt::::new(N1::get_mod());let t2=(m1_3*MInt::::new(N2::get_mod())).inv();let m2=m1*N2::get_mod()as u64;let mut result=vec![0u64;len.min(f.0[0].len())];for(part,((f1,f2),f3))in f.0 .into_iter().zip(f.1).zip(f.2).enumerate(){if f1.is_empty(){continue;}for(value,((c1,c2),c3))in result.iter_mut().zip(Convolve::::inverse_transform(f1,len).into_iter().zip(Convolve::::inverse_transform(f2,len)).zip(Convolve::::inverse_transform(f3,len))){let d1=c1.inner();let d2=((c2-MInt::::from(d1))*t1).inner();let x=MInt::::new(d1)+MInt::::new(d2)*m1_3;let d3=((c3-x)*t2).inner();let limb=(d1 as u64).wrapping_add((d2 as u64).wrapping_mul(m1)).wrapping_add((d3 as u64).wrapping_mul(m2));*value=value.wrapping_add(limb<<(part*bits));}}result}fn multiply(f:&mut Self::F,g:&Self::F){fn multiply(f:&mut[MVec;3],g:&[MVec;3]){assert_eq!(f[0].len(),g[0].len());if f[1].is_empty()||g[1].is_empty(){if f[1].is_empty()&&!g[1].is_empty(){f[1]=f[0].clone();Convolve::::multiply(&mut f[1],&g[1]);}else if!f[1].is_empty(){Convolve::::multiply(&mut f[1],&g[0]);}Convolve::::multiply(&mut f[0],&g[0]);return;}#[cfg(target_arch="x86_64")]if use_block_ntt::(f[0].len()){for part in(1..if f[2].is_empty(){2}else{3}).rev(){let mut sum=f[0].clone();Convolve::::multiply(&mut sum,&g[part]);for left in 1..=part{let mut product=f[left].clone();Convolve::::multiply(&mut product,&g[part-left]);for(value,product)in sum.iter_mut().zip(product){*value+=product;}}f[part]=sum;}Convolve::::multiply(&mut f[0],&g[0]);return;}if f[2].is_empty(){for i in 0..f[0].len(){f[1][i]=f[0][i]*g[1][i]+f[1][i]*g[0][i];f[0][i]*=g[0][i];}return;}for i in 0..f[0].len(){f[2][i]=f[0][i]*g[2][i]+f[1][i]*g[1][i]+f[2][i]*g[0][i];f[1][i]=f[0][i]*g[1][i]+f[1][i]*g[0][i];f[0][i]*=g[0][i];}}multiply(&mut f.0,&g.0);multiply(&mut f.1,&g.1);multiply(&mut f.2,&g.2);}fn square(t:Self::T,len:usize)->Self::T{let mut f=Self::transform(t,len);let g=f.clone();Self::multiply(&mut f,&g);Self::inverse_transform(f,len)}fn convolve(a:Self::T,b:Self::T)->Self::T{let max_len=Self::length(&a).max(Self::length(&b));let min_len=Self::length(&a).min(Self::length(&b));let(balanced,short)=crate::avx_helper!(@dispatch_avx2_fma(300,64),(1536,512));if max_len<=balanced||min_len<=short{let a_wrapping:&[Wrapping]=unsafe{std::slice::from_raw_parts(a.as_ptr().cast(),a.len())};let b_wrapping:&[Wrapping]=unsafe{std::slice::from_raw_parts(b.as_ptr().cast(),b.len())};let mut c=std::mem::ManuallyDrop::new(if max_len<=300||min_len>60{convolve_karatsuba(a_wrapping,b_wrapping)}else{convolve_naive::<_,true>(a_wrapping,b_wrapping)});return unsafe{Vec::from_raw_parts(c.as_mut_ptr().cast(),c.len(),c.capacity())};}let len=(Self::length(&a)+Self::length(&b)).saturating_sub(1);let block_len=if min_len>=1<<20{1<<20}else{(min_len.next_power_of_two()*8).min(1<<21)-min_len+1};if max_len<=block_len{return convolve_u64_fft(a,b);}let mut result=vec![0u64;len];for(i,a)in a.chunks(block_len).enumerate(){for(j,b)in b.chunks(block_len).enumerate(){if a.len().min(b.len())<=60{for(x,&a)in a.iter().enumerate(){for(y,&b)in b.iter().enumerate(){let value=&mut result[(i+j)*block_len+x+y];*value=value.wrapping_add(a.wrapping_mul(b));}}continue;}let product=convolve_u64_fft(a.to_vec(),b.to_vec());for(value,product)in result[(i+j)*block_len..].iter_mut().zip(product){*value=value.wrapping_add(product);}}}result}}fn convolve_u64_fft(a:Vec,b:Vec)->Vec{crate::avx_helper!(@dispatch_avx2_fma return unsafe{convolve_u64_avx2(a,b)},());convolve_u64_fft_scalar(a,b)}fn convolve_u64_fft_scalar(a:Vec,b:Vec)->Vec{fn split(values:&[u64])->[Vec;5]{let mut result=std::array::from_fn(|_|Vec::with_capacity(values.len()));for mut value in values.iter().copied(){for part in&mut result{let digit=((value<<51)as i64)>>51;part.push(digit);value=(value>>13).wrapping_add(u64::from(digit<0));}}result}let len=a.len()+b.len()-1;let transform=|values:&[u64]|{if values.iter().any(|&value|value>u32::MAX as u64){return split(values).map(|part|ConvolveRealFft::transform(part,len));}let[a,b,c,_,_]=split(values);let a=ConvolveRealFft::transform(a,len);let size=a.len();[a,ConvolveRealFft::transform(b,len),ConvolveRealFft::transform(c,len),vec![Zero::zero();size],vec![Zero::zero();size]]};let fa=transform(&a);drop(a);let fb=transform(&b);drop(b);let values:[Vec;5]=std::array::from_fn(|part|{let mut sum=fa[0].clone();ConvolveRealFft::multiply(&mut sum,&fb[part]);for left in 1..=part{let mut product=fa[left].clone();ConvolveRealFft::multiply(&mut product,&fb[part-left]);for(sum,product)in sum.iter_mut().zip(product){*sum+=product;}}ConvolveRealFft::inverse_transform(sum,len)});(0..len).map(|i|{(values[0][i]as u64).wrapping_add((values[1][i]as u64)<<13).wrapping_add((values[2][i]as u64)<<26).wrapping_add((values[3][i]as u64)<<39).wrapping_add((values[4][i]as u64)<<52)}).collect()}pub trait NttReuse:ConvolveSteps{const MULTIPLE:bool=true;#[doc=" Transforms coefficients into the usual NTT frequency order."]fn transform_ntt(t:Self::T,len:usize)->Self::F{Self::transform(t,len)}#[doc=" Inverts a value produced by `transform_ntt`."]fn inverse_transform_ntt(f:Self::F,len:usize)->Self::T{Self::inverse_transform(f,len)}#[doc=" Extends a value produced by `transform_ntt` to twice its length."]#[doc=" If `monic`, the input represents a monic degree-`n` polynomial modulo"]#[doc=" `x^n - 1`, where `n` is the transform length."]fn ntt_doubling(f:Self::F,monic:bool)->Self::F;#[doc=" Extracts the even coefficients of `a(x) * b(-x)` in the usual NTT frequency order."]fn even_mul_normal_neg(f:&Self::F,g:&Self::F)->Self::F;#[doc=" Extracts the odd coefficients of `a(x) * b(-x)` in the usual NTT frequency order."]fn odd_mul_normal_neg(f:&Self::F,g:&Self::F)->Self::F;#[doc=" Multiplies a usual NTT transform by the corresponding prefix of another one."]fn multiply_prefix(f:&mut Self::F,g:&Self::F);#[doc=" Adds the pointwise product of two usual NTT transforms to `sum`."]fn multiply_add(sum:&mut Self::F,f:&Self::F,g:&Self::F);#[doc=" Maximum number of products that can be summed before reconstruction."]#[doc=" Both factors must transform canonical coefficients at the supplied transform's length,"]#[doc=" and each cyclic product must itself be reconstructible."]fn max_product_sum_count(_f:&Self::F)->usize{if Self::MULTIPLE{1}else{usize::MAX}}fn power_projection_step(p_flat:Self::T,q_flat:Self::T,n:usize,py:usize,qy:usize)->(Self::T,Self::T){let base=n*2;let len_p=base*py;let len_q=base*qy;let len=(len_p+len_q-1).max(len_q+len_q-1);let half=len.max(1).next_power_of_two()/2;let p_fft=Self::transform_ntt(p_flat,len);let q_fft=Self::transform_ntt(q_flat,len);let pr_fft=Self::odd_mul_normal_neg(&p_fft,&q_fft);let qr_fft=Self::even_mul_normal_neg(&q_fft,&q_fft);(Self::inverse_transform_ntt(pr_fft,half),Self::inverse_transform_ntt(qr_fft,half))}}thread_local!(static BIT_REVERSE:UnsafeCell>> =const{UnsafeCell::new(vec![])};);implNttReuse for Convolvewhere M:Montgomery32NttModulus{const MULTIPLE:bool=false;fn transform_ntt(mut t:Self::T,len:usize)->Self::F{t.resize_with(len.max(1).next_power_of_two(),Zero::zero);ntt(&mut t);t}fn inverse_transform_ntt(mut f:Self::F,len:usize)->Self::T{intt(&mut f);f.truncate(len);f}fn ntt_doubling(mut f:Self::F,monic:bool)->Self::F{let n=f.len();let k=n.trailing_zeros()as usize;let mut a=Self::inverse_transform_ntt(f.clone(),n);if monic{a[0]-=MInt::::from(2);}let zeta=MInt::::new_unchecked(M::INFO.root[k+1]);let zeta2=zeta*zeta;let mut rot=[MInt::one(),zeta,zeta2,zeta2*zeta];let step=zeta2*zeta2;for a in a.chunks_mut(4){for(a,rot)in a.iter_mut().zip(&mut rot){*a*=*rot;*rot*=step;}}f.extend(Self::transform_ntt(a,n));f}fn even_mul_normal_neg(f:&Self::F,g:&Self::F)->Self::F{assert_eq!(f.len(),g.len());assert!(f.len().is_power_of_two());assert!(f.len()>=2);if std::ptr::eq(f,g){return f.as_chunks::<2>().0 .iter().map(|a|a[0]*a[1]).collect();}let inv2=MInt::::from(2).inv();let n=f.len()/2;(0..n).map(|i|(f[i<<1]*g[i<<1|1]+f[i<<1|1]*g[i<<1])*inv2).collect()}fn odd_mul_normal_neg(f:&Self::F,g:&Self::F)->Self::F{assert_eq!(f.len(),g.len());assert!(f.len().is_power_of_two());assert!(f.len()>=2);let mut inv2=MInt::::from(2).inv();let n=f.len()/2;let k=f.len().trailing_zeros()as usize;let mut h=vec![MInt::::zero();n];let w=MInt::::new_unchecked(M::INFO.inv_root[k]);BIT_REVERSE.with(|br|{let br=unsafe{&mut*br.get()};if br.len()>1]>>1)|((i&1)<>,q_flat:Vec>,n:usize,py:usize,qy:usize)->(Vec>,Vec>){let high_degree=(qy-1)*2;let rows=(py+qy-1).max(high_degree).next_power_of_two();let cols=n*2;let size=rows*cols;let mut p=p_flat;p.resize_with(size,MInt::::zero);ntt_rows(&mut p,cols);ntt_batch(&mut p,cols);let mut q=q_flat;q.resize_with(size,MInt::::zero);ntt_rows(&mut q,cols);let q_high=(rows==high_degree).then(| |q[(qy-1)*cols..qy*cols].to_vec());ntt_batch(&mut q,cols);let half=cols/2;let mut odd_factor=vec![MInt::::zero();half];let mut factor=MInt::::from(2).inv();let k=cols.trailing_zeros()as usize;let w=MInt::::new_unchecked(M::INFO.inv_root[k]);BIT_REVERSE.with(|br|{let br=unsafe{&mut*br.get()};if br.len()>1]>>1)|((i&1)<::zero();rows*half];let mut qr=vec![MInt::::zero();rows*half];for i in 0..pr.len(){pr[i]=(p[i<<1]*q[i<<1|1]-p[i<<1|1]*q[i<<1])*odd_factor[i&(half-1)];qr[i]=q[i<<1]*q[i<<1|1];}intt_batch(&mut pr,half);intt_rows(&mut pr,half);intt_batch(&mut qr,half);intt_rows(&mut qr,half);if let Some(q_high)=q_high{let mut q_high_even=vec![MInt::::zero();half];for i in 0..half{q_high_even[i]=q_high[i<<1]*q_high[i<<1|1];}intt(&mut q_high_even);for(value,high)in qr.iter_mut().zip(&q_high_even){*value-=*high;}qr.extend_from_slice(&q_high_even);}(pr,qr)}}implNttReuse for Convolve<(M,(N1,N2,N3))>where M:MIntConvert+MIntConvert,N1:Montgomery32NttModulus,N2:Montgomery32NttModulus,N3:Montgomery32NttModulus{fn max_product_sum_count(f:&Self::F)->usize{let modulus=>::mod_into()as u128;if modulus==1{return usize::MAX;}let capacity=N1::MOD as u128*N2::MOD as u128*N3::MOD as u128;((capacity-1)/((modulus-1)*(modulus-1))/f.0 .len()as u128).clamp(1,usize::MAX as u128)as usize}fn transform_ntt(t:Self::T,len:usize)->Self::F{let npot=len.max(1).next_power_of_two();let f=convert_crt_input(t,npot);(Convolve::::transform_ntt(f.0,npot),Convolve::::transform_ntt(f.1,npot),Convolve::::transform_ntt(f.2,npot))}fn inverse_transform_ntt(f:Self::F,len:usize)->Self::T{reconstruct_mint_crt((Convolve::::inverse_transform_ntt(f.0,len),Convolve::::inverse_transform_ntt(f.1,len),Convolve::::inverse_transform_ntt(f.2,len)))}fn ntt_doubling(f:Self::F,monic:bool)->Self::F{if monic{let n=f.0 .len();let mut coefficients=Self::inverse_transform_ntt(f,n);coefficients[0]-=MInt::::one();coefficients.push(MInt::::one());Self::transform_ntt(coefficients,n*2)}else{(Convolve::::ntt_doubling(f.0,false),Convolve::::ntt_doubling(f.1,false),Convolve::::ntt_doubling(f.2,false))}}fn even_mul_normal_neg(f:&Self::F,g:&Self::F)->Self::F{fn even_mul_normal_neg_corrected(f:&[MInt],g:&[MInt],m:u32)->Vec>where M:Montgomery32NttModulus{let n=f.len();assert_eq!(f.len(),g.len());assert!(f.len().is_power_of_two());assert!(f.len()>=2);let inv2=MInt::::from(2).inv();let u=MInt::::new(m)*MInt::::from(n as u32);let n=f.len()/2;(0..n).map(|i|{(f[i<<1]*if i==0{g[i<<1|1]+u}else{g[i<<1|1]}+f[i<<1|1]*g[i<<1])*inv2}).collect()}let m=M::mod_into();(even_mul_normal_neg_corrected(&f.0,&g.0,m),even_mul_normal_neg_corrected(&f.1,&g.1,m),even_mul_normal_neg_corrected(&f.2,&g.2,m))}fn odd_mul_normal_neg(f:&Self::F,g:&Self::F)->Self::F{fn odd_mul_normal_neg_corrected(f:&[MInt],g:&[MInt],m:u32)->Vec>where M:Montgomery32NttModulus{assert_eq!(f.len(),g.len());assert!(f.len().is_power_of_two());assert!(f.len()>=2);let mut inv2=MInt::::from(2).inv();let u=MInt::::new(m)*MInt::::from(f.len()as u32);let n=f.len()/2;let k=f.len().trailing_zeros()as usize;let mut h=vec![MInt::::zero();n];let w=MInt::::new_unchecked(M::INFO.inv_root[k]);BIT_REVERSE.with(|br|{let br=unsafe{&mut*br.get()};if br.len()>1]>>1)|((i&1)<::multiply_prefix(&mut f.0,&g.0);Convolve::::multiply_prefix(&mut f.1,&g.1);Convolve::::multiply_prefix(&mut f.2,&g.2);}fn multiply_add(sum:&mut Self::F,f:&Self::F,g:&Self::F){Convolve::::multiply_add(&mut sum.0,&f.0,&g.0);Convolve::::multiply_add(&mut sum.1,&f.1,&g.1);Convolve::::multiply_add(&mut sum.2,&f.2,&g.2);}}} // codesnip-guard: PartialIgnoredOrd #[cfg_attr(any(),rustfmt::skip)]pub use self::partial_ignored_ord::PartialIgnoredOrd;#[cfg_attr(any(),rustfmt::skip)]mod partial_ignored_ord{#[derive(Debug,Default,Clone,Copy)]pub struct PartialIgnoredOrd(pub T,pub U);implEq for PartialIgnoredOrd{}implPartialEq for PartialIgnoredOrdwhere T:PartialEq{fn eq(&self,other:&Self)->bool{self.0 .eq(&other.0)}}implPartialOrd for PartialIgnoredOrdwhere T:PartialOrd{fn partial_cmp(&self,other:&Self)->Option{self.0 .partial_cmp(&other.0)}}implOrd for PartialIgnoredOrdwhere T:Ord{fn cmp(&self,other:&Self)->std::cmp::Ordering{self.0 .cmp(&other.0)}}} // codesnip-guard: _huge_pages #[cfg_attr(any(),rustfmt::skip)]pub use self::huge_pages::advise_huge_pages;#[cfg_attr(any(),rustfmt::skip)]mod huge_pages{#[inline]pub fn advise_huge_pages(_values:&mut Vec){#[cfg(all(target_os="linux",target_arch="x86_64"))]unsafe{unsafe extern "C"{fn madvise(addr:*mut std::ffi::c_void,len:usize,advice:std::ffi::c_int)->std::ffi::c_int;}const PAGE_SIZE:usize=1<<12;const MADV_HUGEPAGE:std::ffi::c_int=14;let ptr=_values.as_mut_ptr().cast::();let offset=ptr.align_offset(PAGE_SIZE);let len=(_values.capacity()*size_of::()).saturating_sub(offset)&!(PAGE_SIZE-1);if len>=1<<20{let _=madvise(ptr.add(offset).cast(),len,MADV_HUGEPAGE);}}}} // codesnip-guard: _partial_fft #[cfg_attr(any(),rustfmt::skip)]pub mod partial_fft{#![doc=" Partial transforms modulo X^(blocks * lanes) - i, leaving lane-sized products."]#![doc=" Buffers contain a nonzero power-of-two number of blocks. SIMD functions require"]#![doc=" their declared CPU features. Real inputs pack their two halves into re and im."]use super::{Complex,One};const ROOTS:[Complex;20]=[Complex{re:f64::from_bits(0x0000000000000000),im:f64::from_bits(0xbff0000000000000)},Complex{re:f64::from_bits(0x3fe6a09e667f3bcd),im:f64::from_bits(0xbfe6a09e667f3bcd)},Complex{re:f64::from_bits(0x3fed906bcf328d46),im:f64::from_bits(0xbfd87de2a6aea963)},Complex{re:f64::from_bits(0x3fef6297cff75cb0),im:f64::from_bits(0xbfc8f8b83c69a60b)},Complex{re:f64::from_bits(0x3fefd88da3d12526),im:f64::from_bits(0xbfb917a6bc29b42c)},Complex{re:f64::from_bits(0x3feff621e3796d7e),im:f64::from_bits(0xbfa91f65f10dd814)},Complex{re:f64::from_bits(0x3feffd886084cd0d),im:f64::from_bits(0xbf992155f7a3667e)},Complex{re:f64::from_bits(0x3fefff62169b92db),im:f64::from_bits(0xbf8921d1fcdec784)},Complex{re:f64::from_bits(0x3fefffd8858e8a92),im:f64::from_bits(0xbf7921f0fe670071)},Complex{re:f64::from_bits(0x3feffff621621d02),im:f64::from_bits(0xbf6921f8becca4ba)},Complex{re:f64::from_bits(0x3feffffd88586ee6),im:f64::from_bits(0xbf5921faaee6472e)},Complex{re:f64::from_bits(0x3fefffff62161a34),im:f64::from_bits(0xbf4921fb2aecb360)},Complex{re:f64::from_bits(0x3fefffffd8858675),im:f64::from_bits(0xbf3921fb49ee4ea6)},Complex{re:f64::from_bits(0x3feffffff621619c),im:f64::from_bits(0xbf2921fb51aeb57c)},Complex{re:f64::from_bits(0x3feffffffd885867),im:f64::from_bits(0xbf1921fb539ecf31)},Complex{re:f64::from_bits(0x3fefffffff62161a),im:f64::from_bits(0xbf0921fb541ad59e)},Complex{re:f64::from_bits(0x3fefffffffd88586),im:f64::from_bits(0xbef921fb5439d73a)},Complex{re:f64::from_bits(0x3feffffffff62162),im:f64::from_bits(0xbee921fb544197a1)},Complex{re:f64::from_bits(0x3feffffffffd8858),im:f64::from_bits(0xbed921fb544387ba)},Complex{re:f64::from_bits(0x3fefffffffff6216),im:f64::from_bits(0xbec921fb544403c1)}];thread_local!{static ROTATE_CACHE:std::cell::RefCell>> =std::cell::RefCell::new(vec![Complex::one()]);}#[doc=" Borrows at least `n` cached roots, where `n` is a nonzero power of two."]#[doc=" The callback must not reenter this cache."]pub fn with_roots(n:usize,f:impl FnOnce(&[Complex])->T)->T{assert!(n.is_power_of_two());ROTATE_CACHE.with_borrow_mut(|cache|{let mut m=cache.len();if m],step:usize,n:usize,k:usize)->Complex{let k=step*k;let w=cache[(k>>2)<<1].conjugate();let w=match k&3{0=>w,1=>Complex::new(-w.re,-w.im),2=>Complex::new(-w.im,w.re),_=>Complex::new(w.im,-w.re),};cache[step*n].conjugate()*w}macro_rules!partial_fft{([$($qualifier:tt)*]$(#[$attr:meta])*;$block:ident,$lanes:literal)=>{$(#[$attr])*pub$($qualifier)*fn fft_soa(a:&mut[$block]){let n=a.len()*$lanes;with_roots(n/2,|cache|{let parity=(n.trailing_zeros()-($lanes as u32).trailing_zeros())&1;for leaf in(0..n).step_by(4*$lanes){let mut level=(n+leaf).trailing_zeros();level-=u32::from(level&1!=n.trailing_zeros()&1);while level>=($lanes as u32).trailing_zeros()+2{let len=1usize<>level;let width=len/(4*$lanes);let start=q*width*4;let(a,rest)=a[start..start+width*4].split_at_mut(width);let(b,rest)=rest.split_at_mut(width);let(c,d)=rest.split_at_mut(width);let w1=eval_twiddle(cache,4,n>>level,q);let w2=w1*w1;let w3=w1*w2;let(w1r,w1i)=(splat(w1.re),splat(w1.im));let(w2r,w2i)=(splat(w2.re),splat(w2.im));let(w3r,w3i)=(splat(w3.re),splat(w3.im));for i in 0..width{let(ar,ai)=load(&a[i]);let(br,bi)=load(&b[i]);let(cr,ci)=load(&c[i]);let(dr,di)=load(&d[i]);let(br,bi)=complex_mul(br,bi,w1r,w1i);let(cr,ci)=complex_mul(cr,ci,w2r,w2i);let(dr,di)=complex_mul(dr,di,w3r,w3i);let acr=add(ar,cr);let aci=add(ai,ci);let bdr=add(br,dr);let bdi=add(bi,di);let acd_r=sub(ar,cr);let acd_i=sub(ai,ci);let bdd_r=sub(br,dr);let bdd_i=sub(bi,di);store(&mut a[i],add(acr,bdr),add(aci,bdi));store(&mut b[i],sub(acr,bdr),sub(aci,bdi));store(&mut c[i],sub(acd_r,bdd_i),add(acd_i,bdd_r));store(&mut d[i],add(acd_r,bdd_i),sub(acd_i,bdd_r));}level-=2;}}if parity!=0{let blocks=n/(2*$lanes);for k in 0..blocks{let w=eval_twiddle(cache,2,blocks,k);let wr=splat(w.re);let wi=splat(w.im);let(ar,ai)=load(&a[k*2]);let(br,bi)=load(&a[k*2+1]);let(br,bi)=complex_mul(br,bi,wr,wi);store(&mut a[k*2],add(ar,br),add(ai,bi));store(&mut a[k*2+1],sub(ar,br),sub(ai,bi));}}});}$(#[$attr])*pub$($qualifier)*fn ifft_soa(a:&mut[$block]){let n=a.len()*$lanes;with_roots(n/2,|cache|{let parity=(n.trailing_zeros()-($lanes as u32).trailing_zeros())&1;if parity!=0{let blocks=n/(2*$lanes);for k in 0..blocks{let w=eval_twiddle(cache,2,blocks,k).conjugate();let wr=splat(w.re);let wi=splat(w.im);let(ar,ai)=load(&a[k*2]);let(br,bi)=load(&a[k*2+1]);store(&mut a[k*2],add(ar,br),add(ai,bi));let(br,bi)=complex_mul(sub(ar,br),sub(ai,bi),wr,wi);store(&mut a[k*2+1],br,bi);}}for leaf in(3*$lanes..n).step_by(4*$lanes){let max_level=(leaf+$lanes-1).trailing_ones();let mut level=($lanes as u32).trailing_zeros()+2+parity;while level<=max_level{let len=1usize<>level;let width=len/(4*$lanes);let start=q*width*4;let(a,rest)=a[start..start+width*4].split_at_mut(width);let(b,rest)=rest.split_at_mut(width);let(c,d)=rest.split_at_mut(width);let w1=eval_twiddle(cache,4,n>>level,q).conjugate();let w2=w1*w1;let w3=w1*w2;let(w1r,w1i)=(splat(w1.re),splat(w1.im));let(w2r,w2i)=(splat(w2.re),splat(w2.im));let(w3r,w3i)=(splat(w3.re),splat(w3.im));for i in 0..width{let(ar,ai)=load(&a[i]);let(br,bi)=load(&b[i]);let(cr,ci)=load(&c[i]);let(dr,di)=load(&d[i]);let abr=add(ar,br);let abi=add(ai,bi);let cdr=add(cr,dr);let cdi=add(ci,di);let abd_r=sub(ar,br);let abd_i=sub(ai,bi);let cdd_r=sub(cr,dr);let cdd_i=sub(ci,di);store(&mut a[i],add(abr,cdr),add(abi,cdi));let(br,bi)=complex_mul(add(abd_r,cdd_i),sub(abd_i,cdd_r),w1r,w1i);store(&mut b[i],br,bi);let(cr,ci)=complex_mul(sub(abr,cdr),sub(abi,cdi),w2r,w2i);store(&mut c[i],cr,ci);let(dr,di)=complex_mul(sub(abd_r,cdd_i),add(abd_i,cdd_r),w3r,w3i);store(&mut d[i],dr,di);}level+=2;}}let scale=splat($lanes as f64/n as f64);for value in a{let(re,im)=load(value);store(value,mul(re,scale),mul(im,scale));}});}};}pub mod scalar{use super::{eval_twiddle,with_roots};#[derive(Clone,Copy,Default)]pub struct Complex4{pub re:[f64;4],pub im:[f64;4]}#[inline]pub fn load(v:&Complex4)->([f64;4],[f64;4]){(v.re,v.im)}#[inline]pub fn store(v:&mut Complex4,re:[f64;4],im:[f64;4]){v.re=re;v.im=im;}#[inline]pub fn splat(x:f64)->[f64;4]{[x;4]}#[inline]pub fn first(a:f64,b:f64)->[f64;4]{[a,b,b,b]}#[inline]pub fn rotate(v:[f64;4])->[f64;4]{[v[3],v[0],v[1],v[2]]}#[inline]pub fn add(a:[f64;4],b:[f64;4])->[f64;4]{std::array::from_fn(|i|a[i]+b[i])}#[inline]pub fn sub(a:[f64;4],b:[f64;4])->[f64;4]{std::array::from_fn(|i|a[i]-b[i])}#[inline]pub fn mul(a:[f64;4],b:[f64;4])->[f64;4]{std::array::from_fn(|i|a[i]*b[i])}#[inline]pub fn complex_mul(ar:[f64;4],ai:[f64;4],br:[f64;4],bi:[f64;4])->([f64;4],[f64;4]){(sub(mul(ar,br),mul(ai,bi)),add(mul(ai,br),mul(ar,bi)))}#[inline]pub fn multiply_accumulate(rr:&mut[f64;4],ri:&mut[f64;4],ar:[f64;4],ai:[f64;4],br:[f64;4],bi:[f64;4]){*rr=add(*rr,mul(ar,br));*rr=sub(*rr,mul(ai,bi));*ri=add(*ri,mul(ai,br));*ri=add(*ri,mul(ar,bi));}partial_fft!([];Complex4,4);}#[cfg(target_arch="x86_64")]pub mod avx2{#![allow(clippy::missing_safety_doc,unsafe_op_in_unsafe_fn)]#[cfg(target_arch="x86_64")]use super::{eval_twiddle,with_roots};#[cfg(target_arch="x86_64")]use std::arch::x86_64::*;#[cfg(target_arch="x86_64")]pub use std::arch::x86_64::{_mm256_add_pd as add,_mm256_mul_pd as mul,_mm256_set1_pd as splat,_mm256_sub_pd as sub};#[derive(Clone,Copy,Default)]#[repr(C,align(32))]#[cfg(target_arch="x86_64")]pub struct Complex4{pub re:[f64;4],pub im:[f64;4]}#[target_feature(enable="avx2,fma")]#[inline]#[cfg(target_arch="x86_64")]pub unsafe fn load(value:&Complex4)->(__m256d,__m256d){(_mm256_load_pd(value.re.as_ptr()),_mm256_load_pd(value.im.as_ptr()))}#[target_feature(enable="avx2,fma")]#[inline]#[cfg(target_arch="x86_64")]pub unsafe fn store(value:&mut Complex4,re:__m256d,im:__m256d){_mm256_store_pd(value.re.as_mut_ptr(),re);_mm256_store_pd(value.im.as_mut_ptr(),im);}#[target_feature(enable="avx2,fma")]#[inline]#[cfg(target_arch="x86_64")]pub unsafe fn complex_mul(ar:__m256d,ai:__m256d,br:__m256d,bi:__m256d)->(__m256d,__m256d){(_mm256_fmsub_pd(ar,br,_mm256_mul_pd(ai,bi)),_mm256_fmadd_pd(ai,br,_mm256_mul_pd(ar,bi)))}#[target_feature(enable="avx2,fma")]#[inline]#[cfg(target_arch="x86_64")]pub unsafe fn multiply_accumulate(rr:&mut __m256d,ri:&mut __m256d,ar:__m256d,ai:__m256d,br:__m256d,bi:__m256d){*rr=_mm256_fmadd_pd(ar,br,*rr);*rr=_mm256_fnmadd_pd(ai,bi,*rr);*ri=_mm256_fmadd_pd(ai,br,*ri);*ri=_mm256_fmadd_pd(ar,bi,*ri);}#[target_feature(enable="avx2,fma")]#[inline]#[cfg(target_arch="x86_64")]pub unsafe fn first(a:f64,b:f64)->__m256d{_mm256_setr_pd(a,b,b,b)}#[target_feature(enable="avx2,fma")]#[inline]#[cfg(target_arch="x86_64")]pub unsafe fn rotate(v:__m256d)->__m256d{_mm256_permute4x64_pd::<0x93>(v)}#[cfg(target_arch="x86_64")]partial_fft!([unsafe]#[target_feature(enable="avx2,fma")];Complex4,4);}#[cfg(target_arch="x86_64")]pub mod avx512{#![allow(clippy::missing_safety_doc,unsafe_op_in_unsafe_fn)]#[cfg(target_arch="x86_64")]use super::{eval_twiddle,with_roots};#[cfg(target_arch="x86_64")]use std::arch::x86_64::*;#[cfg(target_arch="x86_64")]pub use std::arch::x86_64::{_mm512_add_pd as add,_mm512_mul_pd as mul,_mm512_set1_pd as splat,_mm512_sub_pd as sub};#[derive(Clone,Copy,Default)]#[repr(C,align(64))]#[cfg(target_arch="x86_64")]pub struct Complex8{pub re:[f64;8],pub im:[f64;8]}#[target_feature(enable="avx2,avx512f,fma")]#[inline]#[cfg(target_arch="x86_64")]pub unsafe fn load(value:&Complex8)->(__m512d,__m512d){(_mm512_load_pd(value.re.as_ptr()),_mm512_load_pd(value.im.as_ptr()))}#[target_feature(enable="avx2,avx512f,fma")]#[inline]#[cfg(target_arch="x86_64")]pub unsafe fn store(value:&mut Complex8,re:__m512d,im:__m512d){_mm512_store_pd(value.re.as_mut_ptr(),re);_mm512_store_pd(value.im.as_mut_ptr(),im);}#[target_feature(enable="avx2,avx512f,fma")]#[inline]#[cfg(target_arch="x86_64")]pub unsafe fn complex_mul(ar:__m512d,ai:__m512d,br:__m512d,bi:__m512d)->(__m512d,__m512d){(_mm512_fmsub_pd(ar,br,_mm512_mul_pd(ai,bi)),_mm512_fmadd_pd(ai,br,_mm512_mul_pd(ar,bi)))}#[target_feature(enable="avx2,avx512f,fma")]#[inline]#[cfg(target_arch="x86_64")]pub unsafe fn multiply_accumulate(rr:&mut __m512d,ri:&mut __m512d,ar:__m512d,ai:__m512d,br:__m512d,bi:__m512d){*rr=_mm512_fmadd_pd(ar,br,*rr);*rr=_mm512_fnmadd_pd(ai,bi,*rr);*ri=_mm512_fmadd_pd(ai,br,*ri);*ri=_mm512_fmadd_pd(ar,bi,*ri);}#[target_feature(enable="avx2,avx512f,fma")]#[inline]#[cfg(target_arch="x86_64")]pub unsafe fn first(a:f64,b:f64)->__m512d{_mm512_setr_pd(a,b,b,b,b,b,b,b)}#[target_feature(enable="avx2,avx512f,fma")]#[inline]#[cfg(target_arch="x86_64")]pub unsafe fn rotate(v:__m512d)->__m512d{_mm512_permutexvar_pd(_mm512_setr_epi64(7,0,1,2,3,4,5,6),v)}#[cfg(target_arch="x86_64")]partial_fft!([unsafe]#[target_feature(enable="avx2,avx512f,fma")];Complex8,8);}} // codesnip-guard: array #[cfg_attr(any(),rustfmt::skip)]mod array{#[macro_export]macro_rules!array{[@inner$data:ident=[$init:expr;$len:expr]]=>{{use::std::mem::{ManuallyDrop,MaybeUninit};let mut$data:[MaybeUninit<_>;$len]=unsafe{MaybeUninit::uninit().assume_init()};$init;#[repr(C)]union __Transmuter{src:ManuallyDrop<[MaybeUninit;N]>,dst:ManuallyDrop<[T;N]>,}ManuallyDrop::into_inner(unsafe{__Transmuter{src:ManuallyDrop::new($data)}.dst})}};[||$e:expr;$len:expr]=>{$crate::array![@inner data=[data.iter_mut().for_each(|item|*item=MaybeUninit::new($e));$len]]};[|$i:pat_param|$e:expr;$len:expr]=>{$crate::array![@inner data=[data.iter_mut().enumerate().for_each(|($i,item)|*item=MaybeUninit::new($e));$len]]};[$e:expr;$len:expr]=>{{let e=$e;$crate::array![||Clone::clone(&e);$len]}};}} // 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: bounded #[cfg_attr(any(),rustfmt::skip)]pub use self::bounded::Bounded;#[cfg_attr(any(),rustfmt::skip)]mod bounded{#[doc=" Trait for max/min bounds"]pub trait Bounded:Sized+PartialOrd{fn maximum()->Self;fn minimum()->Self;fn is_maximum(&self)->bool{self==&Self::maximum()}fn is_minimum(&self)->bool{self==&Self::minimum()}fn set_maximum(&mut self){*self=Self::maximum()}fn set_minimum(&mut self){*self=Self::minimum()}}macro_rules!impl_bounded_num{($($t:ident)*)=>{$(impl Bounded for$t{fn maximum()->Self{$t::MAX}fn minimum()->Self{$t::MIN}})*};}impl_bounded_num!(u8 u16 u32 u64 u128 usize i8 i16 i32 i64 i128 isize f32 f64);macro_rules!impl_bounded_tuple{(@impl$($T:ident)*)=>{impl<$($T:Bounded),*>Bounded for($($T,)*){fn maximum()->Self{($(<$T as Bounded>::maximum(),)*)}fn minimum()->Self{($(<$T as Bounded>::minimum(),)*)}}};(@inner$($T:ident)*,)=>{impl_bounded_tuple!(@impl$($T)*);};(@inner$($T:ident)*,$U:ident$($Rest:ident)*)=>{impl_bounded_tuple!(@impl$($T)*);impl_bounded_tuple!(@inner$($T)*$U,$($Rest)*);};($T:ident$($Rest:ident)*)=>{impl_bounded_tuple!(@inner$T,$($Rest)*);};}impl_bounded_tuple!(A B C D E F G H I J);impl Bounded for(){fn maximum()->Self{}fn minimum()->Self{}}impl Bounded for bool{fn maximum()->Self{true}fn minimum()->Self{false}}implBounded for Optionwhere T:Bounded{fn maximum()->Self{Some(::maximum())}fn minimum()->Self{None}}implBounded for std::cmp::Reversewhere T:Bounded{fn maximum()->Self{std::cmp::Reverse(::minimum())}fn minimum()->Self{std::cmp::Reverse(::maximum())}}} // codesnip-guard: factorial #[cfg_attr(any(),rustfmt::skip)]pub use self::factorial::MemorizedFactorial;#[cfg_attr(any(),rustfmt::skip)]mod factorial{use super::{MInt,MIntConvert,One,Zero};#[derive(Clone,Debug)]pub struct MemorizedFactorialwhere M:MIntConvert{pub fact:Vec>,pub inv_fact:Vec>}implMemorizedFactorialwhere M:MIntConvert{pub fn new(max_n:usize)->Self{let mut fact=vec![MInt::one();max_n+1];let mut inv_fact=vec![MInt::one();max_n+1];for i in 2..=max_n{fact[i]=fact[i-1]*MInt::from(i);}inv_fact[max_n]=fact[max_n].inv();for i in(3..=max_n).rev(){inv_fact[i-1]=inv_fact[i]*MInt::from(i);}Self{fact,inv_fact}}pub fn combination(&self,n:usize,r:usize)->MInt{debug_assert!(nMInt{debug_assert!(nMInt{debug_assert!(n+rMInt{debug_assert!(n0);self.inv_fact[n]*self.fact[n-1]}}} // codesnip-guard: fast_fourier_transform #[cfg_attr(any(),rustfmt::skip)]pub use self::fast_fourier_transform::ConvolveRealFft;#[cfg_attr(any(),rustfmt::skip)]mod fast_fourier_transform{use super::partial_fft::with_roots;use super::{Complex,ConvolveSteps,Zero};#[cfg(target_arch="x86_64")]use super::{SimdBackend,advise_huge_pages,simd_backend};pub enum ConvolveRealFft{}#[cfg(target_arch="x86_64")]pub mod simd{#![allow(clippy::missing_safety_doc,unsafe_op_in_unsafe_fn)]#[cfg(target_arch="x86_64")]use super::{advise_huge_pages,with_roots};#[cfg(target_arch="x86_64")]use std::arch::x86_64::*;#[cfg(target_arch="x86_64")]use super::super::partial_fft::{avx2::{Complex4,complex_mul,fft_soa,ifft_soa,load,multiply_accumulate,store},eval_twiddle};#[target_feature(enable="avx2,fma")]#[cfg(target_arch="x86_64")]unsafe fn dot_one_soa(a:&mut[Complex4],b:&[Complex4]){let n=a.len()*4;with_roots(n/2,|cache|{for i in 0..a.len(){let(mut br,mut bi)=load(&b[i]);let mut rr=_mm256_setzero_pd();let mut ri=_mm256_setzero_pd();let w=eval_twiddle(cache,1,a.len(),i);let wr=_mm256_setr_pd(w.re,1.0,1.0,1.0);let wi=_mm256_setr_pd(w.im,0.0,0.0,0.0);for lane in 0..4{let ar=_mm256_set1_pd(a[i].re[lane]);let ai=_mm256_set1_pd(a[i].im[lane]);multiply_accumulate(&mut rr,&mut ri,ar,ai,br,bi);if lane!=3{br=_mm256_permute4x64_pd::<0x93>(br);bi=_mm256_permute4x64_pd::<0x93>(bi);(br,bi)=complex_mul(br,bi,wr,wi);}}store(&mut a[i],rr,ri);}});}#[inline]#[cfg(target_arch="x86_64")]fn pack_f64(values:impl Iterator,n:usize)->Vec{let mut result=Vec::with_capacity(n/4);advise_huge_pages(&mut result);result.resize(n/4,Complex4::default());for(i,value)in values.enumerate(){if i>2].re[i&3]=value;}else{result[(i-n)>>2].im[i&3]=value;}}result}#[target_feature(enable="avx2,fma")]#[cfg(target_arch="x86_64")]pub unsafe fn convolve_f64_avx2(a:impl ExactSizeIterator,b:impl ExactSizeIterator,range:std::ops::Range)->Vec{let n=(range.end.next_power_of_two()/2).max(4);let mut fa=pack_f64(a,n);let mut fb=pack_f64(b,n);fft_soa(&mut fa);fft_soa(&mut fb);dot_one_soa(&mut fa,&fb);drop(fb);ifft_soa(&mut fa);range.map(|i|{if i>2].re[i&3]}else{fa[(i-n)>>2].im[i&3]}}).collect()}#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]pub unsafe fn convolve_i64_avx2(a:Vec,b:Vec,len:usize)->Vec{super::convolve_i64_naive(a,b,len)}#[target_feature(enable="avx512f,avx512dq,avx512cd,avx512bw,avx512vl")]#[cfg(target_arch="x86_64")]pub unsafe fn convolve_i64_avx512(a:Vec,b:Vec,len:usize)->Vec{super::convolve_i64_naive(a,b,len)}}fn bit_reverse(f:&mut[T]){let mut ip=vec![0u32];let mut k=f.len();let mut m=1;while 2*m)){const BLOCK:usize=256;let sign=if inverse{1.0}else{-1.0};let step=Complex::primitive_nth_root_of_unity(sign*n as f64);for start in(1..n/4).step_by(BLOCK){let mut w=Complex::polar(1.0,sign*std::f64::consts::TAU*start as f64/n as f64);for k in start..(start+BLOCK).min(n/4){f(k,w);w*=step;}}}pub fn transform_real(t:impl IntoIterator,len:usize)->Vec>{let n=len.max(4).next_power_of_two();let mut f=vec![Complex::zero();n/2];for(i,t)in t.into_iter().enumerate(){if i&1==0{f[i/2].re=t;}else{f[i/2].im=t;}}fft(&mut f);bit_reverse(&mut f);f[0]=Complex::new(f[0].re+f[0].im,f[0].re-f[0].im);f[n/4]=f[n/4].conjugate();real_twiddles(n,false,|k,wk|{let c=wk.conjugate().transpose()+1.;let d=c*(f[k]-f[n/2-k].conjugate())*0.5;f[k]-=d;f[n/2-k]+=d.conjugate();});f}pub fn inverse_transform_real(mut f:Vec>,len:usize)->Vec{let n=len.max(4).next_power_of_two();assert_eq!(f.len(),n/2);f[0]=Complex::new((f[0].re+f[0].im)*0.5,(f[0].re-f[0].im)*0.5);f[n/4]=f[n/4].conjugate();real_twiddles(n,true,|k,wk|{let c=wk.transpose().conjugate()+1.;let d=c*(f[k]-f[n/2-k].conjugate())*0.5;f[k]-=d;f[n/2-k]+=d.conjugate();});bit_reverse(&mut f);ifft(&mut f);let inv=1./(n/2)as f64;(0..len).map(|i|inv*if i&1==0{f[i/2].re}else{f[i/2].im}).collect()}#[inline(always)]fn convolve_i64_naive(a:Vec,b:Vec,len:usize)->Vec{let(a,b)=if a.len();type F=Vec>;fn length(t:&Self::T)->usize{t.len()}fn transform(t:Self::T,len:usize)->Self::F{transform_real(t.into_iter().map(|t|t as f64),len)}fn inverse_transform(f:Self::F,len:usize)->Self::T{inverse_transform_real(f,len).into_iter().map(|value|value.round()as i64).collect()}fn convolve(a:Self::T,b:Self::T)->Self::T{let len=(a.len()+b.len()).saturating_sub(1);if(a.len().min(b.len())<=32||{let size=len.next_power_of_two();let log=size.ilog2();let limit=crate::avx_helper!(@dispatch simd_backend,SimdBackend;3*log+16,2*log+16,4*log+16);2*a.len()as u128*b.len()as u128<=limit as u128*size as u128})&&a.iter().map(|x|x.unsigned_abs()).max().unwrap_or(0)as u128*b.iter().map(|x|x.unsigned_abs()).max().unwrap_or(0)as u128<=(1u128<<53)/a.len().min(b.len()).max(1)as u128{return crate::avx_helper!(@dispatch simd_backend,SimdBackend;unsafe{if a.len().max(b.len())<32{simd::convolve_i64_avx2(a,b,len)}else{simd::convolve_i64_avx512(a,b,len)}},unsafe{simd::convolve_i64_avx2(a,b,len)},convolve_i64_naive(a,b,len));}if!a.is_empty()&&!b.is_empty(){crate::avx_helper!(@dispatch_avx2_fma return unsafe{simd::convolve_f64_avx2(a.into_iter().map(|x|x as f64),b.into_iter().map(|x|x as f64),0..len,).into_iter().map(|x|x.round()as i64).collect()},());}let mut a=Self::transform(a,len);let b=Self::transform(b,len);Self::multiply(&mut a,&b);Self::inverse_transform(a,len)}fn multiply(f:&mut Self::F,g:&Self::F){assert_eq!(f.len(),g.len());f[0].re*=g[0].re;f[0].im*=g[0].im;for(f,g)in f.iter_mut().zip(g.iter()).skip(1){*f*=*g;}}}fn middle_product_f64_scalar(a:impl ExactSizeIterator,b:impl ExactSizeIterator)->Vec{let a_len=a.len();let b_len=b.len();let len=a_len+b_len-1;let mut a=transform_real(a,len);let b=transform_real(b,len);ConvolveRealFft::multiply(&mut a,&b);inverse_transform_real(a,len)[b_len-1..a_len].to_vec()}impl ConvolveRealFft{#[doc=" Returns coefficients `b.len() - 1..a.len()` of the convolution of `a` and `b`."]#[doc=" Panics unless `0 < b.len() <= a.len()`."]pub fn middle_product_f64(a:impl ExactSizeIterator,b:impl ExactSizeIterator)->Vec{assert!(0{{let a=$a;let cache=$cache;let n=a.len();if$inverse{let mut v=1;if n.trailing_zeros()&1==1{for(a,w)in a.as_chunks_mut::<2>().0 .iter_mut().zip(cache){let y=(a[0]-a[1])*w.conjugate();a[0]+=a[1];a[1]=y;}v=2;}while v=2{let l=v/2;for(q,block)in a.chunks_exact_mut(l*4).enumerate(){let(a,rest)=block.split_at_mut(l);let(b,rest)=rest.split_at_mut(l);let(c,d)=rest.split_at_mut(l);let w0=cache[q];let w1=cache[q<<1];let w3=w0*w1;for(((a,b),c),d)in a.iter_mut().zip(b).zip(c).zip(d){let bv=*b*w1;let cv=*c*w0;let dv=*d*w3;let ac0=*a+cv;let ac1=*a-cv;let bd0=bv+dv;let bd1=bv-dv;let bd1=Complex::new(bd1.im,-bd1.re);*a=ac0+bd0;*b=ac0-bd0;*c=ac1+bd1;*d=ac1-bd1;}}v>>=2;}if v==1{for(a,w)in a.as_chunks_mut::<2>().0 .iter_mut().zip(cache){let y=a[1]**w;a[1]=a[0]-y;a[0]+=y;}}}}};}pub fn fft(a:&mut[Complex]){fft_dispatch::(a);}pub fn ifft(a:&mut[Complex]){fft_dispatch::(a);}fn fft_dispatch(a:&mut[Complex]){with_roots(a.len()/2,|cache|{#[cfg(target_arch="x86_64")]if a.len()>=16{match simd_backend(){SimdBackend::Avx512=>{return unsafe{fft_avx512::(a,cache)};}SimdBackend::Avx2=>return unsafe{fft_avx2::(a,cache)},SimdBackend::Scalar=>{}}}fft_kernel!(a,cache,INVERSE);});}#[cfg(target_arch="x86_64")]#[target_feature(enable="avx2")]unsafe fn fft_avx2(a:&mut[Complex],cache:&[Complex]){fft_kernel!(a,cache,INVERSE);}#[cfg(target_arch="x86_64")]#[target_feature(enable="avx512f,avx512dq,avx512cd,avx512bw,avx512vl")]unsafe fn fft_avx512(a:&mut[Complex],cache:&[Complex]){fft_kernel!(a,cache,INVERSE);}} // codesnip-guard: float #[cfg_attr(any(),rustfmt::skip)]pub use self::float::{Float,Float32,Float64};#[cfg_attr(any(),rustfmt::skip)]mod float{use super::{Bounded,One,Scan,ScanSource,Zero};use std::{cmp::Ordering,fmt::Display,num::FpCategory,ops::{Add,Div,Mul,Neg,Rem,Sub},str::FromStr};pub trait Float:Copy+Default+Display+FromStr+PartialEq+PartialOrd+Zero+One+Bounded+Add+Sub+Mul+Div+Neg+Rem{type Bits;fn floor(self)->Self;fn ceil(self)->Self;fn round(self)->Self;fn trunc(self)->Self;fn fract(self)->Self;fn abs(self)->Self;fn signum(self)->Self;fn copysign(self,sign:Self)->Self;fn mul_add(self,a:Self,b:Self)->Self;fn div_euclid(self,rhs:Self)->Self;fn rem_euclid(self,rhs:Self)->Self;fn powi(self,n:i32)->Self;fn powf(self,n:Self)->Self;fn sqrt(self)->Self;fn exp(self)->Self;fn exp2(self)->Self;fn ln(self)->Self;fn log(self,base:Self)->Self;fn log2(self)->Self;fn log10(self)->Self;fn cbrt(self)->Self;fn hypot(self,other:Self)->Self;fn sin(self)->Self;fn cos(self)->Self;fn tan(self)->Self;fn asin(self)->Self;fn acos(self)->Self;fn atan(self)->Self;fn atan2(self,other:Self)->Self;fn sin_cos(self)->(Self,Self);fn exp_m1(self)->Self;fn ln_1p(self)->Self;fn sinh(self)->Self;fn cosh(self)->Self;fn tanh(self)->Self;fn asinh(self)->Self;fn acosh(self)->Self;fn atanh(self)->Self;fn is_nan(self)->bool;fn is_infinite(self)->bool;fn is_finite(self)->bool;fn is_normal(self)->bool;fn classify(self)->FpCategory;fn is_sign_positive(self)->bool;fn is_sign_negative(self)->bool;fn recip(self)->Self;fn to_degrees(self)->Self;fn to_radians(self)->Self;fn max(self,other:Self)->Self;fn min(self,other:Self)->Self;fn midpoint(self,other:Self)->Self;fn to_bits(self)->Self::Bits;fn from_bits(v:Self::Bits)->Self;fn total_cmp(&self,other:&Self)->Ordering;fn clamp(self,min:Self,max:Self)->Self;const RADIX:u32;const MANTISSA_DIGITS:u32;const DIGITS:u32;const EPSILON:Self;const MIN:Self;const MIN_POSITIVE:Self;const MAX:Self;const MIN_EXP:i32;const MAX_EXP:i32;const MIN_10_EXP:i32;const MAX_10_EXP:i32;const NAN:Self;const INFINITY:Self;const NEG_INFINITY:Self;const PI:Self;const TAU:Self;const FRAC_PI_2:Self;const FRAC_PI_3:Self;const FRAC_PI_4:Self;const FRAC_PI_6:Self;const FRAC_PI_8:Self;const FRAC_1_PI:Self;const FRAC_2_PI:Self;const FRAC_2_SQRT_PI:Self;const SQRT_2:Self;const FRAC_1_SQRT_2:Self;const E:Self;const LOG2_E:Self;const LOG10_E:Self;const LN_2:Self;const LN_10:Self;}macro_rules!impl_primitive_float{($({$t:ident$i:ident$u:ident$e:expr})*)=>{$(impl Float for$t{type Bits=$u;fn floor(self)->Self{$t::floor(self)}fn ceil(self)->Self{$t::ceil(self)}fn round(self)->Self{$t::round(self)}fn trunc(self)->Self{$t::trunc(self)}fn fract(self)->Self{$t::fract(self)}fn abs(self)->Self{$t::abs(self)}fn signum(self)->Self{$t::signum(self)}fn copysign(self,sign:Self)->Self{$t::copysign(self,sign)}fn mul_add(self,a:Self,b:Self)->Self{$t::mul_add(self,a,b)}fn div_euclid(self,rhs:Self)->Self{$t::div_euclid(self,rhs)}fn rem_euclid(self,rhs:Self)->Self{$t::rem_euclid(self,rhs)}fn powi(self,n:i32)->Self{$t::powi(self,n)}fn powf(self,n:Self)->Self{$t::powf(self,n)}fn sqrt(self)->Self{$t::sqrt(self)}fn exp(self)->Self{$t::exp(self)}fn exp2(self)->Self{$t::exp2(self)}fn ln(self)->Self{$t::ln(self)}fn log(self,base:Self)->Self{$t::log(self,base)}fn log2(self)->Self{$t::log2(self)}fn log10(self)->Self{$t::log10(self)}fn cbrt(self)->Self{$t::cbrt(self)}fn hypot(self,other:Self)->Self{$t::hypot(self,other)}fn sin(self)->Self{$t::sin(self)}fn cos(self)->Self{$t::cos(self)}fn tan(self)->Self{$t::tan(self)}fn asin(self)->Self{$t::asin(self)}fn acos(self)->Self{$t::acos(self)}fn atan(self)->Self{$t::atan(self)}fn atan2(self,other:Self)->Self{$t::atan2(self,other)}fn sin_cos(self)->(Self,Self){$t::sin_cos(self)}fn exp_m1(self)->Self{$t::exp_m1(self)}fn ln_1p(self)->Self{$t::ln_1p(self)}fn sinh(self)->Self{$t::sinh(self)}fn cosh(self)->Self{$t::cosh(self)}fn tanh(self)->Self{$t::tanh(self)}fn asinh(self)->Self{$t::asinh(self)}fn acosh(self)->Self{$t::acosh(self)}fn atanh(self)->Self{$t::atanh(self)}fn is_nan(self)->bool{$t::is_nan(self)}fn is_infinite(self)->bool{$t::is_infinite(self)}fn is_finite(self)->bool{$t::is_finite(self)}fn is_normal(self)->bool{$t::is_normal(self)}fn classify(self)->std::num::FpCategory{$t::classify(self)}fn is_sign_positive(self)->bool{$t::is_sign_positive(self)}fn is_sign_negative(self)->bool{$t::is_sign_negative(self)}fn recip(self)->Self{$t::recip(self)}fn to_degrees(self)->Self{$t::to_degrees(self)}fn to_radians(self)->Self{$t::to_radians(self)}fn max(self,other:Self)->Self{$t::max(self,other)}fn min(self,other:Self)->Self{$t::min(self,other)}fn midpoint(self,other:Self)->Self{$t::midpoint(self,other)}fn to_bits(self)->Self::Bits{$t::to_bits(self)}fn from_bits(v:Self::Bits)->Self{$t::from_bits(v)}fn total_cmp(&self,other:&Self)->Ordering{$t::total_cmp(self,other)}fn clamp(self,min:Self,max:Self)->Self{$t::clamp(self,min,max)}const RADIX:u32=$t::RADIX;const MANTISSA_DIGITS:u32=$t::MANTISSA_DIGITS;const DIGITS:u32=$t::DIGITS;const EPSILON:Self=$t::EPSILON;const MIN:Self=$t::MIN;const MIN_POSITIVE:Self=$t::MIN_POSITIVE;const MAX:Self=$t::MAX;const MIN_EXP:i32=$t::MIN_EXP;const MAX_EXP:i32=$t::MAX_EXP;const MIN_10_EXP:i32=$t::MIN_10_EXP;const MAX_10_EXP:i32=$t::MAX_10_EXP;const NAN:Self=$t::NAN;const INFINITY:Self=$t::INFINITY;const NEG_INFINITY:Self=$t::NEG_INFINITY;const PI:Self=std::$t::consts::PI;const TAU:Self=std::$t::consts::PI*2.0;const FRAC_PI_2:Self=std::$t::consts::FRAC_PI_2;const FRAC_PI_3:Self=std::$t::consts::FRAC_PI_3;const FRAC_PI_4:Self=std::$t::consts::FRAC_PI_4;const FRAC_PI_6:Self=std::$t::consts::FRAC_PI_6;const FRAC_PI_8:Self=std::$t::consts::FRAC_PI_8;const FRAC_1_PI:Self=std::$t::consts::FRAC_1_PI;const FRAC_2_PI:Self=std::$t::consts::FRAC_2_PI;const FRAC_2_SQRT_PI:Self=std::$t::consts::FRAC_2_SQRT_PI;const SQRT_2:Self=std::$t::consts::SQRT_2;const FRAC_1_SQRT_2:Self=std::$t::consts::FRAC_1_SQRT_2;const E:Self=std::$t::consts::E;const LOG2_E:Self=std::$t::consts::LOG2_E;const LOG10_E:Self=std::$t::consts::LOG10_E;const LN_2:Self=std::$t::consts::LN_2;const LN_10:Self=std::$t::consts::LN_10;})*};}impl_primitive_float!({f32 i32 u32 31}{f64 i64 u64 63});macro_rules!impl_ord_float{($({$t:ident$n:ident})*)=>{$(#[derive(Debug,Copy,Clone,Default)]#[repr(transparent)]pub struct$n(pub$t);impl std::fmt::Display for$n{fn fmt(&self,f:&mut std::fmt::Formatter<'_>)->std::fmt::Result{<$t as std::fmt::Display>::fmt(&self.0,f)}}impl std::str::FromStr for$n{type Err=std::num::ParseFloatError;fn from_str(s:&str)->Result{<$t as std::str::FromStr>::from_str(s).map(Self)}}impl From<$t>for$n{fn from(x:$t)->Self{Self(x)}}impl Zero for$n{fn zero()->Self{Self(<$t as Zero>::zero())}fn is_zero(&self)->bool{self.0==0.0}}impl One for$n{fn one()->Self{Self(<$t as One>::one())}}impl Bounded for$n{fn maximum()->Self{Self(<$t as Bounded>::maximum())}fn minimum()->Self{Self(<$t as Bounded>::minimum())}}impl Add for$n{type Output=Self;fn add(self,rhs:Self)->Self::Output{Self(<$t as Add>::add(self.0,rhs.0))}}impl Sub for$n{type Output=Self;fn sub(self,rhs:Self)->Self::Output{Self(<$t as Sub>::sub(self.0,rhs.0))}}impl Mul for$n{type Output=Self;fn mul(self,rhs:Self)->Self::Output{Self(<$t as Mul>::mul(self.0,rhs.0))}}impl Div for$n{type Output=Self;fn div(self,rhs:Self)->Self::Output{Self(<$t as Div>::div(self.0,rhs.0))}}impl Neg for$n{type Output=Self;fn neg(self)->Self::Output{Self(<$t as Neg>::neg(self.0))}}impl Rem for$n{type Output=Self;fn rem(self,rhs:Self)->Self::Output{Self(<$t as Rem>::rem(self.0,rhs.0))}}impl PartialEq for$n{fn eq(&self,other:&Self)->bool{self.total_cmp(other).is_eq()}}impl Eq for$n{}impl PartialOrd for$n{fn partial_cmp(&self,other:&Self)->Option{Some(self.cmp(other))}}impl Ord for$n{fn cmp(&self,other:&Self)->Ordering{self.total_cmp(other)}}impl Scan for$n{type Output=Self;fn scan(iter:&mut I)->Option{<$t as Scan>::scan(iter).map(Self)}}impl Float for$n{type Bits=<$t as Float>::Bits;fn floor(self)->Self{Self(<$t as Float>::floor(self.0))}fn ceil(self)->Self{Self(<$t as Float>::ceil(self.0))}fn round(self)->Self{Self(<$t as Float>::round(self.0))}fn trunc(self)->Self{Self(<$t as Float>::trunc(self.0))}fn fract(self)->Self{Self(<$t as Float>::fract(self.0))}fn abs(self)->Self{Self(<$t as Float>::abs(self.0))}fn signum(self)->Self{Self(<$t as Float>::signum(self.0))}fn copysign(self,sign:Self)->Self{Self(<$t as Float>::copysign(self.0,sign.0))}fn mul_add(self,a:Self,b:Self)->Self{Self(<$t as Float>::mul_add(self.0,a.0,b.0))}fn div_euclid(self,rhs:Self)->Self{Self(<$t as Float>::div_euclid(self.0,rhs.0))}fn rem_euclid(self,rhs:Self)->Self{Self(<$t as Float>::rem_euclid(self.0,rhs.0))}fn powi(self,n:i32)->Self{Self(<$t as Float>::powi(self.0,n))}fn powf(self,n:Self)->Self{Self(<$t as Float>::powf(self.0,n.0))}fn sqrt(self)->Self{Self(<$t as Float>::sqrt(self.0))}fn exp(self)->Self{Self(<$t as Float>::exp(self.0))}fn exp2(self)->Self{Self(<$t as Float>::exp2(self.0))}fn ln(self)->Self{Self(<$t as Float>::ln(self.0))}fn log(self,base:Self)->Self{Self(<$t as Float>::log(self.0,base.0))}fn log2(self)->Self{Self(<$t as Float>::log2(self.0))}fn log10(self)->Self{Self(<$t as Float>::log10(self.0))}fn cbrt(self)->Self{Self(<$t as Float>::cbrt(self.0))}fn hypot(self,other:Self)->Self{Self(<$t as Float>::hypot(self.0,other.0))}fn sin(self)->Self{Self(<$t as Float>::sin(self.0))}fn cos(self)->Self{Self(<$t as Float>::cos(self.0))}fn tan(self)->Self{Self(<$t as Float>::tan(self.0))}fn asin(self)->Self{Self(<$t as Float>::asin(self.0))}fn acos(self)->Self{Self(<$t as Float>::acos(self.0))}fn atan(self)->Self{Self(<$t as Float>::atan(self.0))}fn atan2(self,other:Self)->Self{Self(<$t as Float>::atan2(self.0,other.0))}fn sin_cos(self)->(Self,Self){let(sin,cos)=<$t as Float>::sin_cos(self.0);(Self(sin),Self(cos))}fn exp_m1(self)->Self{Self(<$t as Float>::exp_m1(self.0))}fn ln_1p(self)->Self{Self(<$t as Float>::ln_1p(self.0))}fn sinh(self)->Self{Self(<$t as Float>::sinh(self.0))}fn cosh(self)->Self{Self(<$t as Float>::cosh(self.0))}fn tanh(self)->Self{Self(<$t as Float>::tanh(self.0))}fn asinh(self)->Self{Self(<$t as Float>::asinh(self.0))}fn acosh(self)->Self{Self(<$t as Float>::acosh(self.0))}fn atanh(self)->Self{Self(<$t as Float>::atanh(self.0))}fn is_nan(self)->bool{<$t as Float>::is_nan(self.0)}fn is_infinite(self)->bool{<$t as Float>::is_infinite(self.0)}fn is_finite(self)->bool{<$t as Float>::is_finite(self.0)}fn is_normal(self)->bool{<$t as Float>::is_normal(self.0)}fn classify(self)->std::num::FpCategory{<$t as Float>::classify(self.0)}fn is_sign_positive(self)->bool{<$t as Float>::is_sign_positive(self.0)}fn is_sign_negative(self)->bool{<$t as Float>::is_sign_negative(self.0)}fn recip(self)->Self{Self(<$t as Float>::recip(self.0))}fn to_degrees(self)->Self{Self(<$t as Float>::to_degrees(self.0))}fn to_radians(self)->Self{Self(<$t as Float>::to_radians(self.0))}fn max(self,other:Self)->Self{Self(<$t as Float>::max(self.0,other.0))}fn min(self,other:Self)->Self{Self(<$t as Float>::min(self.0,other.0))}fn midpoint(self,other:Self)->Self{Self(<$t as Float>::midpoint(self.0,other.0))}fn to_bits(self)->Self::Bits{<$t as Float>::to_bits(self.0)}fn from_bits(v:Self::Bits)->Self{Self(<$t as Float>::from_bits(v))}fn total_cmp(&self,other:&Self)->Ordering{<$t as Float>::total_cmp(&self.0,&other.0)}fn clamp(self,min:Self,max:Self)->Self{Self(<$t as Float>::clamp(self.0,min.0,max.0))}const RADIX:u32=<$t as Float>::RADIX;const MANTISSA_DIGITS:u32=<$t as Float>::MANTISSA_DIGITS;const DIGITS:u32=<$t as Float>::DIGITS;const EPSILON:Self=Self(<$t as Float>::EPSILON);const MIN:Self=Self(<$t as Float>::MIN);const MIN_POSITIVE:Self=Self(<$t as Float>::MIN_POSITIVE);const MAX:Self=Self(<$t as Float>::MAX);const MIN_EXP:i32=<$t as Float>::MIN_EXP;const MAX_EXP:i32=<$t as Float>::MAX_EXP;const MIN_10_EXP:i32=<$t as Float>::MIN_10_EXP;const MAX_10_EXP:i32=<$t as Float>::MAX_10_EXP;const NAN:Self=Self(<$t as Float>::NAN);const INFINITY:Self=Self(<$t as Float>::INFINITY);const NEG_INFINITY:Self=Self(<$t as Float>::NEG_INFINITY);const PI:Self=Self(<$t as Float>::PI);const TAU:Self=Self(<$t as Float>::TAU);const FRAC_PI_2:Self=Self(<$t as Float>::FRAC_PI_2);const FRAC_PI_3:Self=Self(<$t as Float>::FRAC_PI_3);const FRAC_PI_4:Self=Self(<$t as Float>::FRAC_PI_4);const FRAC_PI_6:Self=Self(<$t as Float>::FRAC_PI_6);const FRAC_PI_8:Self=Self(<$t as Float>::FRAC_PI_8);const FRAC_1_PI:Self=Self(<$t as Float>::FRAC_1_PI);const FRAC_2_PI:Self=Self(<$t as Float>::FRAC_2_PI);const FRAC_2_SQRT_PI:Self=Self(<$t as Float>::FRAC_2_SQRT_PI);const SQRT_2:Self=Self(<$t as Float>::SQRT_2);const FRAC_1_SQRT_2:Self=Self(<$t as Float>::FRAC_1_SQRT_2);const E:Self=Self(<$t as Float>::E);const LOG2_E:Self=Self(<$t as Float>::LOG2_E);const LOG10_E:Self=Self(<$t as Float>::LOG10_E);const LN_2:Self=Self(<$t as Float>::LN_2);const LN_10:Self=Self(<$t as Float>::LN_10);})*};}impl_ord_float!({f32 Float32}{f64 Float64});} // codesnip-guard: gcd #[cfg_attr(any(),rustfmt::skip)]#[doc=" binary gcd"]pub fn gcd(mut a:u64,mut b:u64)->u64{if a==0{return b;}if b==0{return a;}let u=a.trailing_zeros();let v=b.trailing_zeros();a>>=u;b>>=v;let k=std::cmp::min(u,v);while a!=b{if a>=a.trailing_zeros();}a<MIntwhere M:MIntConvert{pub fn sqrt(self)->Option{fn jacobi(mut x:u32)->i8 where M:MIntConvert{let mut s=1i8;let mut m=M::mod_into();while m>1{x%=m;if x==0{return 0;}let k=x.trailing_zeros();if k%2==1&&(m+2)&4!=0{s=-s;}x>>=k;if x&m&2!=0{s=-s;}std::mem::swap(&mut x,&mut m);}s}if M::mod_into()==2{return Some(self);}let j=jacobi::(u32::from(self));match j.cmp(&0){Ordering::Less=>{return None;}Ordering::Equal=>{return Some(Self::zero());}Ordering::Greater=>{}}let mut r=1;let(mut f0,d)=loop{r^=r<<5;r^=r>>17;r^=r<<11;let b=Self::from(r);let d=b*b-self;if jacobi::(u32::from(d))==-1{break(b,d);}};let(mut f1,mut g0,mut g1,mut e)=(Self::one(),Self::one(),Self::zero(),M::mod_into().div_ceil(2));while e>0{if e%2==1{let t=g0*f0+d*g1*f1;g1=g0*f1+g1*f0;g0=t;}let t=f0*f0+d*f1*f1;f1=Self::from(2)*f0*f1;f0=t;e/=2;}if u32::from(g0)>M::mod_into()-u32::from(g0){g0=-g0;}Some(g0)}}} // codesnip-guard: montgomery #[cfg_attr(any(),rustfmt::skip)]pub mod montgomery{use super::*;implMIntBase for M where M:MontgomeryReduction32{type Inner=u32;fn get_mod()->Self::Inner{::MOD}fn mod_zero()->Self::Inner{0}fn mod_one()->Self::Inner{Self::N1}fn mod_add(x:Self::Inner,y:Self::Inner)->Self::Inner{let z=x+y;let m=Self::get_mod();if z>=m{z-m}else{z}}fn mod_sub(x:Self::Inner,y:Self::Inner)->Self::Inner{if xSelf::Inner{Self::reduce(x as u64*y as u64)}fn mod_div(x:Self::Inner,y:Self::Inner)->Self::Inner{Self::mod_mul(x,Self::mod_inv(y))}fn mod_neg(x:Self::Inner)->Self::Inner{if x==0{0}else{Self::get_mod()-x}}fn mod_inv(x:Self::Inner)->Self::Inner{let p=Self::get_mod()as i32;let(mut a,mut b)=(x as i32,p);let(mut u,mut x)=(1,0);while a!=0{let k=b/a;x-=k*u;b-=k*a;std::mem::swap(&mut x,&mut u);std::mem::swap(&mut b,&mut a);}Self::reduce((if x<0{x+p}else{x})as u64*Self::N3 as u64)}fn mod_inner(x:Self::Inner)->Self::Inner{let y=((x.wrapping_mul(Self::R.wrapping_neg())as u64*Self::MOD as u64)>>32)as u32;if y==0{0}else{Self::MOD-y}}}implMIntConvertfor M where M:MontgomeryReduction32{fn from(x:u32)->Self::Inner{Self::reduce(x as u64*Self::N2 as u64)}fn into(x:Self::Inner)->u32{Self::mod_inner(x)}fn mod_into()->u32{::get_mod()}}implMIntConvertfor M where M:MontgomeryReduction32{fn from(x:u64)->Self::Inner{Self::reduce(x%Self::get_mod()as u64*Self::N2 as u64)}fn into(x:Self::Inner)->u64{Self::mod_inner(x)as u64}fn mod_into()->u64{::get_mod()as u64}}implMIntConvertfor M where M:MontgomeryReduction32{fn from(x:usize)->Self::Inner{Self::reduce(x as u64%Self::get_mod()as u64*Self::N2 as u64)}fn into(x:Self::Inner)->usize{Self::mod_inner(x)as usize}fn mod_into()->usize{::get_mod()as usize}}implMIntConvertfor M where M:MontgomeryReduction32{fn from(x:i32)->Self::Inner{let x=x%::get_mod()as i32;let x=if x<0{(x+::get_mod()as i32)as u64}else{x as u64};Self::reduce(x*Self::N2 as u64)}fn into(x:Self::Inner)->i32{Self::mod_inner(x)as i32}fn mod_into()->i32{::get_mod()as i32}}implMIntConvertfor M where M:MontgomeryReduction32{fn from(x:i64)->Self::Inner{let x=x%::get_mod()as i64;let x=if x<0{(x+::get_mod()as i64)as u64}else{x as u64};Self::reduce(x*Self::N2 as u64)}fn into(x:Self::Inner)->i64{Self::mod_inner(x)as i64}fn mod_into()->i64{::get_mod()as i64}}implMIntConvertfor M where M:MontgomeryReduction32{fn from(x:isize)->Self::Inner{let x=x%::get_mod()as isize;let x=if x<0{(x+::get_mod()as isize)as u64}else{x as u64};Self::reduce(x*Self::N2 as u64)}fn into(x:Self::Inner)->isize{Self::mod_inner(x)as isize}fn mod_into()->isize{::get_mod()as isize}}#[doc=" m is prime, n = 2^32"]pub trait MontgomeryReduction32{#[doc=" m"]const MOD:u32;#[doc=" (-m)^{-1} mod n"]const R:u32={let m=Self::MOD;let mut r=0;let mut t=0;let mut i=0;while i<32{if t%2==0{t+=m;r+=1<u32{let m:u32=Self::MOD;let r=Self::R;let mut x=((x+r.wrapping_mul(x as u32)as u64*m as u64)>>32)as u32;if x>=m{x-=m;}x}}macro_rules!define_montgomery_reduction_32{($([$name:ident,$m:expr,$mint_name:ident$(,)?]),*$(,)?)=>{$(pub enum$name{}impl MontgomeryReduction32 for$name{const MOD:u32=$m;}pub type$mint_name=MInt<$name>;)*};}define_montgomery_reduction_32!([Modulo167772161,167_772_161,MInt167772161],[Modulo469762049,469_762_049,MInt469762049],[Modulo754974721,754_974_721,MInt754974721],[Modulo998244353,998_244_353,MInt998244353],);} // codesnip-guard: montgomery_simd #[cfg_attr(any(),rustfmt::skip)]#[cfg(target_arch="x86_64")]pub mod montgomery_simd{#![allow(unsafe_op_in_unsafe_fn)]#[cfg(target_arch="x86_64")]use std::arch::x86_64::*;#[inline]#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]pub unsafe fn montgomery_mul_256(a:__m256i,b:__m256i,r_vec:__m256i,mod_vec:__m256i)->__m256i{let a13=_mm256_bsrli_epi128::<4>(a);let b13=_mm256_bsrli_epi128::<4>(b);let t02=_mm256_mul_epu32(a,b);let t13=_mm256_mul_epu32(a13,b13);let m02=mul_epu32_256(t02,r_vec);let m13=mul_epu32_256(t13,r_vec);let u02=_mm256_add_epi64(t02,_mm256_mul_epu32(m02,mod_vec));let u13=_mm256_add_epi64(t13,_mm256_mul_epu32(m13,mod_vec));_mm256_or_si256(_mm256_bsrli_epi128::<4>(u02),u13)}#[inline]#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]pub unsafe fn montgomery_mul_256_fixed(a:__m256i,b:__m256i,b_r:__m256i,mod_vec:__m256i)->__m256i{let a13=_mm256_bsrli_epi128::<4>(a);let t02=_mm256_mul_epu32(a,b);let t13=_mm256_mul_epu32(a13,b);let m02=mul_epu32_256(a,b_r);let m13=mul_epu32_256(a13,b_r);let u02=_mm256_add_epi64(t02,_mm256_mul_epu32(m02,mod_vec));let u13=_mm256_add_epi64(t13,_mm256_mul_epu32(m13,mod_vec));_mm256_or_si256(_mm256_bsrli_epi128::<4>(u02),u13)}#[inline]#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]pub unsafe fn add_mod_256(a:__m256i,b:__m256i,mod_vec:__m256i)->__m256i{let sum=_mm256_add_epi32(a,b);_mm256_min_epu32(sum,_mm256_sub_epi32(sum,mod_vec))}#[inline]#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]pub unsafe fn sub_mod_256(a:__m256i,b:__m256i,mod_vec:__m256i)->__m256i{let diff=_mm256_sub_epi32(_mm256_add_epi32(a,mod_vec),b);_mm256_min_epu32(diff,_mm256_sub_epi32(diff,mod_vec))}#[inline]#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]pub unsafe fn montgomery_mul_256_canon(a:__m256i,b:__m256i,r_vec:__m256i,mod_vec:__m256i)->__m256i{let x=montgomery_mul_256(a,b,r_vec,mod_vec);_mm256_min_epu32(x,_mm256_sub_epi32(x,mod_vec))}#[inline]#[target_feature(enable="avx512f,avx512dq,avx512cd,avx512bw,avx512vl")]#[cfg(target_arch="x86_64")]pub unsafe fn montgomery_mul_512(a:__m512i,b:__m512i,r_vec:__m512i,mod_vec:__m512i)->__m512i{let a13=_mm512_srli_epi64::<32>(a);let b13=_mm512_srli_epi64::<32>(b);let t02=_mm512_mul_epu32(a,b);let t13=_mm512_mul_epu32(a13,b13);let m02=_mm512_mul_epu32(t02,r_vec);let m13=_mm512_mul_epu32(t13,r_vec);let u02=_mm512_add_epi64(t02,_mm512_mul_epu32(m02,mod_vec));let u13=_mm512_add_epi64(t13,_mm512_mul_epu32(m13,mod_vec));_mm512_or_si512(_mm512_srli_epi64::<32>(u02),u13)}#[inline]#[target_feature(enable="avx512f,avx512dq,avx512cd,avx512bw,avx512vl")]#[cfg(target_arch="x86_64")]pub unsafe fn add_mod_512(a:__m512i,b:__m512i,mod_vec:__m512i)->__m512i{let sum=_mm512_add_epi32(a,b);_mm512_min_epu32(sum,_mm512_sub_epi32(sum,mod_vec))}#[inline]#[target_feature(enable="avx512f,avx512dq,avx512cd,avx512bw,avx512vl")]#[cfg(target_arch="x86_64")]pub unsafe fn sub_mod_512(a:__m512i,b:__m512i,mod_vec:__m512i)->__m512i{let diff=_mm512_sub_epi32(_mm512_add_epi32(a,mod_vec),b);_mm512_min_epu32(diff,_mm512_sub_epi32(diff,mod_vec))}#[inline]#[target_feature(enable="avx512f,avx512dq,avx512cd,avx512bw,avx512vl")]#[cfg(target_arch="x86_64")]pub unsafe fn montgomery_mul_512_canon(a:__m512i,b:__m512i,r_vec:__m512i,mod_vec:__m512i)->__m512i{let x=montgomery_mul_512(a,b,r_vec,mod_vec);_mm512_min_epu32(x,_mm512_sub_epi32(x,mod_vec))}#[inline]#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]pub unsafe fn mul_epu32_256(a:__m256i,b:__m256i)->__m256i{let product;std::arch::asm!("vpmuludq {product}, {a}, {b}",a=in(ymm_reg)a,b=in(ymm_reg)b,product=lateout(ymm_reg)product,options(pure,nomem,nostack,preserves_flags),);product}#[inline]#[target_feature(enable="avx512f")]#[cfg(target_arch="x86_64")]pub unsafe fn montgomery_mul_512_fixed(a:__m512i,b:__m512i,b_r:__m512i,mod_vec:__m512i)->__m512i{let a13=_mm512_srli_epi64::<32>(a);let t02=_mm512_mul_epu32(a,b);let t13=_mm512_mul_epu32(a13,b);let m02=mul_epu32_512(a,b_r);let m13=mul_epu32_512(a13,b_r);let u02=_mm512_add_epi64(t02,_mm512_mul_epu32(m02,mod_vec));let u13=_mm512_add_epi64(t13,_mm512_mul_epu32(m13,mod_vec));_mm512_or_si512(_mm512_srli_epi64::<32>(u02),u13)}#[inline]#[target_feature(enable="avx512f")]#[cfg(target_arch="x86_64")]pub unsafe fn mul_epu32_512(a:__m512i,b:__m512i)->__m512i{let product;std::arch::asm!("vpmuludq {product}, {a}, {b}",a=in(zmm_reg)a,b=in(zmm_reg)b,product=lateout(zmm_reg)product,options(pure,nomem,nostack,preserves_flags),);product}} // codesnip-guard: scanner #[cfg_attr(any(),rustfmt::skip)]pub use self::scanner::*;#[cfg_attr(any(),rustfmt::skip)]mod scanner{use std::{iter::{FromIterator,from_fn,repeat_with},marker::PhantomData};pub fn read_stdin_all()->String{use std::io::Read as _;let mut s=String::new();std::io::stdin().read_to_string(&mut s).expect("io error");s}pub fn read_stdin_all_unchecked()->String{use std::io::Read as _;let mut buf=Vec::new();std::io::stdin().read_to_end(&mut buf).expect("io error");unsafe{String::from_utf8_unchecked(buf)}}pub fn read_all(mut reader:impl std::io::Read)->String{let mut s=String::new();reader.read_to_string(&mut s).expect("io error");s}pub fn read_all_unchecked(mut reader:impl std::io::Read)->String{let mut buf=Vec::new();reader.read_to_end(&mut buf).expect("io error");unsafe{String::from_utf8_unchecked(buf)}}pub fn read_stdin_line()->String{let mut s=String::new();std::io::stdin().read_line(&mut s).expect("io error");s}pub trait Scan:Sized{type Output;fn scan(source:&mut S)->Option;fn scan_vec(source:&mut S,size:usize,capacity:usize)->Vec{let mut values=Vec::with_capacity(capacity.max(size));values.extend((0..size).map(|_|source.scan::()));values}}pub trait MarkedScan:Sized{type Output;fn mscan(self,source:&mut S)->Option;}pub trait ScanSource:Sized{#[doc=" Reads a token under the source's input requirements."]#[doc=" Checked sources return `None` at EOF; unchecked sources require available input."]fn next_token(&mut self)->Option<&str>;#[doc=" Skips separators without consuming the next token."]#[inline]fn skip_whitespace(&mut self){}#[inline]fn read_u8(&mut self)->Option{self.next_token()?.parse().ok()}#[inline]fn read_u16(&mut self)->Option{self.next_token()?.parse().ok()}#[inline]fn read_u32(&mut self)->Option{self.next_token()?.parse().ok()}fn read_u32_vec(&mut self,size:usize,capacity:usize)->Vec{let mut values=Vec::with_capacity(capacity.max(size));values.extend((0..size).map(|_|self.scan::()));values}#[inline]fn read_u64(&mut self)->Option{self.next_token()?.parse().ok()}#[inline]fn read_u128(&mut self)->Option{self.next_token()?.parse().ok()}#[inline]fn read_usize(&mut self)->Option{self.next_token()?.parse().ok()}#[inline]fn read_i8(&mut self)->Option{self.next_token()?.parse().ok()}#[inline]fn read_i16(&mut self)->Option{self.next_token()?.parse().ok()}#[inline]fn read_i32(&mut self)->Option{self.next_token()?.parse().ok()}#[inline]fn read_i64(&mut self)->Option{self.next_token()?.parse().ok()}#[inline]fn read_i128(&mut self)->Option{self.next_token()?.parse().ok()}#[inline]fn read_isize(&mut self)->Option{self.next_token()?.parse().ok()}#[doc=" Panics if reading fails."]#[inline]fn scan(&mut self)->T::Output{T::scan(self).expect("scan error")}#[doc=" Panics if reading fails."]#[inline]fn mscan(&mut self,marker:T)->T::Output{marker.mscan(self).expect("scan error")}fn scan_vec(&mut self,size:usize,capacity:usize)->Vec{if size==0{self.skip_whitespace();}T::scan_vec(self,size,capacity)}#[inline]fn iter(&mut self)->ScannerIter<'_,Self,T>{ScannerIter{inner:self,_marker:PhantomData}}}#[derive(Clone,Debug)]pub struct Scanner<'a,I:Iterator=std::str::SplitAsciiWhitespace<'a>>{iter:I}impl<'a>Scanner<'a>{pub fn new(s:&'a str)->Self{Self{iter:s.split_ascii_whitespace()}}}impl<'a,I:Iterator>Scanner<'a,I>{pub fn new_from_iter(iter:I)->Self{Self{iter}}}impl<'a,I:Iterator>ScanSource for Scanner<'a,I>{fn next_token(&mut self)->Option<&str>{self.iter.next()}}macro_rules!impl_scan{($($t:ty)*)=>{$(impl Scan for$t{type Output=Self;fn scan(iter:&mut I)->Option{iter.next_token()?.parse::<$t>().ok()}})*};}impl_scan!(char f32 f64 String);macro_rules!impl_integer_scan{($($t:ty=>$read:ident),*$(,)?)=>{$(impl Scan for$t{type Output=Self;#[inline]fn scan(source:&mut S)->Option{source.$read()}})*};}impl_integer_scan!(u8=>read_u8,u16=>read_u16,u64=>read_u64,u128=>read_u128,usize=>read_usize,i8=>read_i8,i16=>read_i16,i32=>read_i32,i64=>read_i64,i128=>read_i128,isize=>read_isize,);impl Scan for u32{type Output=Self;#[inline]fn scan(source:&mut S)->Option{source.read_u32()}fn scan_vec(source:&mut S,size:usize,capacity:usize)->Vec{source.read_u32_vec(size,capacity)}}macro_rules!impl_scan_tuple{(@impl$($T:ident)*)=>{impl<$($T:Scan),*>Scan for($($T,)*){type Output=($(<$T as Scan>::Output,)*);fn scan(_iter:&mut It)->Option{Some(($(<$T as Scan>::scan(_iter)?,)*))}}};(@inner$($T:ident)*,)=>{impl_scan_tuple!(@impl$($T)*);};(@inner$($T:ident)*,$U:ident$($Rest:ident)*)=>{impl_scan_tuple!(@impl$($T)*);impl_scan_tuple!(@inner$($T)*$U,$($Rest)*);};($($T:ident)*)=>{impl_scan_tuple!(@inner,$($T)*);};}impl_scan_tuple!(A B C D E F G H I J K);pub struct ScannerIter<'a,S,T>{inner:&'a mut S,_marker:PhantomDataT>}implIterator for ScannerIter<'_,S,T>{type Item=T::Output;#[inline]fn next(&mut self)->Option{T::scan(self.inner)}}#[doc=" scan a value with Scanner"]#[doc=""]#[doc=" - `scan_value!(scanner, ELEMENT)`"]#[doc=""]#[doc=" ELEMENT :="]#[doc=" - `$ty`: Scan"]#[doc=" - `&str`: borrowed token; the source stays borrowed until its last use"]#[doc=" - `@$expr`: MarkedScan"]#[doc=" - `$ty = $expr`: MarkedScan"]#[doc=" - `[ELEMENT; $expr]`: vector"]#[doc=" - `[$ty; $len; capacity $capacity]`: vector with reserved capacity"]#[doc=" - `[ELEMENT; const $expr]`: array"]#[doc=" - `[ELEMENT]`: iterator"]#[doc=" - `[ELEMENT; iter $expr]`: iterator of the specified length"]#[doc=" - `($(ELEMENT)*,)`: tuple"]#[macro_export]macro_rules!scan_value{(@repeat$scanner:expr,[$($t:tt)*]$len:expr)=>{{$crate::scan_value!(@iter$scanner,[$($t)*]$len).collect::>()}};(@repeat$scanner:expr,[$($t:tt)*])=>{{::std::iter::repeat_with(||$crate::scan_value!(@inner$scanner,[]$($t)*))}};(@iter$scanner:expr,[$($t:tt)*]$len:expr)=>{{let size=$len;if size==0{$scanner.skip_whitespace();}$crate::scan_value!(@repeat$scanner,[$($t)*]).take(size)}};(@array$scanner:expr,[$($t:tt)*]$len:expr)=>{{if$len==0{$scanner.skip_whitespace();}$crate::array![||$crate::scan_value!(@inner$scanner,[]$($t)*);$len]}};(@tuple$scanner:expr,[$([$($args:tt)*])*])=>{($($($args)*,)*)};(@sparen$scanner:expr,[]@$e:expr;$($t:tt)*)=>{$crate::scan_value!(@sparen$scanner,[@$e]$($t)*)};(@sparen$scanner:expr,[]($($tt:tt)*);$($t:tt)*)=>{$crate::scan_value!(@sparen$scanner,[($($tt)*)]$($t)*)};(@sparen$scanner:expr,[][$($tt:tt)*];$($t:tt)*)=>{$crate::scan_value!(@sparen$scanner,[[$($tt)*]]$($t)*)};(@sparen$scanner:expr,[]$ty:ty=$e:expr;$($t:tt)*)=>{$crate::scan_value!(@sparen$scanner,[$ty=$e]$($t)*)};(@sparen$scanner:expr,[]$ty:ty;$($t:tt)*)=>{$crate::scan_value!(@sparen$scanner,[$ty]$($t)*)};(@sparen$scanner:expr,[]$($args:tt)*)=>{$crate::scan_value!(@repeat$scanner,[$($args)*])};(@sparen$scanner:expr,[$($args:tt)+]const$len:expr)=>{$crate::scan_value!(@array$scanner,[$($args)+]$len)};(@sparen$scanner:expr,[$($args:tt)+]iter$len:expr)=>{$crate::scan_value!(@iter$scanner,[$($args)+]$len)};(@sparen$scanner:expr,[$($args:tt)+]$len:expr;capacity$capacity:expr)=>{{let size=$len;let capacity=$capacity;$scanner.scan_vec::<$($args)+>(size,capacity)}};(@sparen$scanner:expr,[$($args:tt)+]$len:expr)=>{$crate::scan_value!(@repeat$scanner,[$($args)+]$len)};(@$tag:ident$scanner:expr,[[$($args:tt)*]])=>{$($args)*};(@$tag:ident$scanner:expr,[$($args:tt)*]@$e:expr$(,$($t:tt)*)?)=>{$crate::scan_value!(@$tag$scanner,[$($args)*[$scanner.mscan($e)]]$(,$($t)*)?)};(@$tag:ident$scanner:expr,[$($args:tt)*]($($tuple:tt)*)$($t:tt)*)=>{$crate::scan_value!(@$tag$scanner,[$($args)*[$crate::scan_value!(@tuple$scanner,[]$($tuple)*)]]$($t)*)};(@$tag:ident$scanner:expr,[$($args:tt)*][$($tt:tt)*]$($t:tt)*)=>{$crate::scan_value!(@$tag$scanner,[$($args)*[$crate::scan_value!(@sparen$scanner,[]$($tt)*)]]$($t)*)};(@$tag:ident$scanner:expr,[$($args:tt)*]&str$(,$($t:tt)*)?)=>{$crate::scan_value!(@$tag$scanner,[$($args)*[$scanner.next_token().expect("scan error")]]$(,$($t)*)?)};(@$tag:ident$scanner:expr,[$($args:tt)*]$ty:ty=$e:expr$(,$($t:tt)*)?)=>{$crate::scan_value!(@$tag$scanner,[$($args)*[{let _tmp:$ty=$scanner.mscan($e);_tmp}]]$(,$($t)*)?)};(@$tag:ident$scanner:expr,[$($args:tt)*]$ty:ty$(,$($t:tt)*)?)=>{$crate::scan_value!(@$tag$scanner,[$($args)*[$scanner.scan::<$ty>()]]$(,$($t)*)?)};(@$tag:ident$scanner:expr,[$($args:tt)*],$($t:tt)*)=>{$crate::scan_value!(@$tag$scanner,[$($args)*]$($t)*)};(@$tag:ident$scanner:expr,[$($args:tt)*])=>{::std::compile_error!(::std::stringify!($($args)*))};(src=$src:expr,$($t:tt)*)=>{{let mut __scanner=Scanner::new($src);$crate::scan_value!(@inner __scanner,[]$($t)*)}};(iter=$iter:expr,$($t:tt)*)=>{{let mut __scanner=Scanner::new_from_iter($iter);$crate::scan_value!(@inner __scanner,[]$($t)*)}};($scanner:expr,$($t:tt)*)=>{$crate::scan_value!(@inner$scanner,[]$($t)*)}}#[doc=" scan and bind values with Scanner"]#[doc=""]#[doc=" - `scan!(scanner, $($pat $(: ELEMENT)?),*)`"]#[macro_export]macro_rules!scan{(@assert$p:pat)=>{};(@assert$($p:tt)*)=>{::std::compile_error!(::std::concat!("expected pattern, found `",::std::stringify!($($p)*),"`"));};(@pat$scanner:expr,[][])=>{};(@pat$scanner:expr,[][],$($t:tt)*)=>{$crate::scan!(@pat$scanner,[][]$($t)*)};(@pat$scanner:expr,[$($p:tt)*][]$x:ident$($t:tt)*)=>{$crate::scan!(@pat$scanner,[$($p)*$x][]$($t)*)};(@pat$scanner:expr,[$($p:tt)*][]::$($t:tt)*)=>{$crate::scan!(@pat$scanner,[$($p)*::][]$($t)*)};(@pat$scanner:expr,[$($p:tt)*][]&$($t:tt)*)=>{$crate::scan!(@pat$scanner,[$($p)*&][]$($t)*)};(@pat$scanner:expr,[$($p:tt)*][]($($x:tt)*)$($t:tt)*)=>{$crate::scan!(@pat$scanner,[$($p)*($($x)*)][]$($t)*)};(@pat$scanner:expr,[$($p:tt)*][][$($x:tt)*]$($t:tt)*)=>{$crate::scan!(@pat$scanner,[$($p)*[$($x)*]][]$($t)*)};(@pat$scanner:expr,[$($p:tt)*][]{$($x:tt)*}$($t:tt)*)=>{$crate::scan!(@pat$scanner,[$($p)*{$($x)*}][]$($t)*)};(@pat$scanner:expr,[$($p:tt)*][]:$($t:tt)*)=>{$crate::scan!(@ty$scanner,[$($p)*][]$($t)*)};(@pat$scanner:expr,[$($p:tt)*][]$($t:tt)*)=>{$crate::scan!(@let$scanner,[$($p)*][usize]$($t)*)};(@ty$scanner:expr,[$($p:tt)*][$($tt:tt)*]@$e:expr$(,$($t:tt)*)?)=>{$crate::scan!(@let$scanner,[$($p)*][$($tt)*@$e]$(,$($t)*)?)};(@ty$scanner:expr,[$($p:tt)*][$($tt:tt)*]($($x:tt)*)$($t:tt)*)=>{$crate::scan!(@let$scanner,[$($p)*][$($tt)*($($x)*)]$($t)*)};(@ty$scanner:expr,[$($p:tt)*][$($tt:tt)*][$($x:tt)*]$($t:tt)*)=>{$crate::scan!(@let$scanner,[$($p)*][$($tt)*[$($x)*]]$($t)*)};(@ty$scanner:expr,[$($p:tt)*][$($tt:tt)*]&str$(,$($t:tt)*)?)=>{$crate::scan!(@let$scanner,[$($p)*][$($tt)*&str]$(,$($t)*)?)};(@ty$scanner:expr,[$($p:tt)*][$($tt:tt)*]$ty:ty=$e:expr$(,$($t:tt)*)?)=>{$crate::scan!(@let$scanner,[$($p)*][$($tt)*$ty=$e]$(,$($t)*)?)};(@ty$scanner:expr,[$($p:tt)*][$($tt:tt)*]$ty:ty$(,$($t:tt)*)?)=>{$crate::scan!(@let$scanner,[$($p)*][$($tt)*$ty]$(,$($t)*)?)};(@let$scanner:expr,[$($p:tt)*][$($tt:tt)*]$($t:tt)*)=>{$crate::scan!{@assert$($p)*}let$($p)* =$crate::scan_value!($scanner,$($tt)*);$crate::scan!(@pat$scanner,[][]$($t)*)};(src=$src:expr,$($t:tt)*)=>{let mut __scanner=Scanner::new($src);$crate::scan!(@pat __scanner,[][]$($t)*)};(iter=$iter:expr,$($t:tt)*)=>{let mut __scanner=Scanner::new_from_iter($iter);$crate::scan!(@pat __scanner,[][]$($t)*)};($scanner:expr,$($t:tt)*)=>{$crate::scan!(@pat$scanner,[][]$($t)*)}}#[doc=" define enum scan rules"]#[doc=""]#[doc=" # Example"]#[doc=" ```rust"]#[doc=" # use competitive::{define_enum_scan, tools::{CharsWithBase, Scan, ScanSource, Scanner, Usize1}};"]#[doc=" define_enum_scan! {"]#[doc=" enum Query: u8 {"]#[doc=" 0 => Noop,"]#[doc=" 1 => Args { i: Usize1, s: char },"]#[doc=" 9 => Complex { n: usize, c: [(usize, Vec = CharsWithBase('a')); n] },"]#[doc=" }"]#[doc=" }"]#[doc=" ```"]#[macro_export]macro_rules!define_enum_scan{(@field_ty@repeat[$($t:tt)*]$($len:expr)?)=>{Vec<$crate::define_enum_scan!(@field_ty$($t)*)>};(@field_ty@array[$($t:tt)*]$len:expr)=>{[$crate::define_enum_scan!(@field_ty$($t)*);$len]};(@field_ty@tuple[$([$($args:tt)*])*])=>{($($($args)*,)*)};(@field_ty@sparen[]($($tt:tt)*);$($t:tt)*)=>{$crate::define_enum_scan!(@field_ty@sparen[($($tt)*)]$($t)*)};(@field_ty@sparen[][$($tt:tt)*];$($t:tt)*)=>{$crate::define_enum_scan!(@field_ty@sparen[[$($tt)*]]$($t)*)};(@field_ty@sparen[]$ty:ty=$e:expr;$($t:tt)*)=>{$crate::define_enum_scan!(@field_ty@sparen[$ty=$e]$($t)*)};(@field_ty@sparen[]$ty:ty;$($t:tt)*)=>{$crate::define_enum_scan!(@field_ty@sparen[$ty]$($t)*)};(@field_ty@sparen[]$($args:tt)*)=>{$crate::define_enum_scan!(@field_ty@repeat[$($args)*])};(@field_ty@sparen[$($args:tt)+]const$len:expr)=>{$crate::define_enum_scan!(@field_ty@array[$($args)+]$len)};(@field_ty@sparen[$($args:tt)+]$len:expr)=>{$crate::define_enum_scan!(@field_ty@repeat[$($args)+]$len)};(@field_ty@$tag:ident[$($args:tt)*]($($tuple:tt)*)$($t:tt)*)=>{$crate::define_enum_scan!(@field_ty@$tag[$($args)*[$crate::define_enum_scan!(@field_ty@tuple[]$($tuple)*)]]$($t)*)};(@field_ty@$tag:ident[$($args:tt)*][$($tt:tt)*]$($t:tt)*)=>{$crate::define_enum_scan!(@field_ty@$tag[$($args)*[$crate::define_enum_scan!(@field_ty@sparen[]$($tt)*)]]$($t)*)};(@field_ty@$tag:ident[$($args:tt)*]$ty:ty=$e:expr$(,$($t:tt)*)?)=>{$crate::define_enum_scan!(@field_ty@$tag[$($args)*[$ty]]$(,$($t)*)?)};(@field_ty@$tag:ident[$($args:tt)*]$ty:ty$(,$($t:tt)*)?)=>{$crate::define_enum_scan!(@field_ty@$tag[$($args)*[<$ty as Scan>::Output]]$(,$($t)*)?)};(@field_ty@$tag:ident[$($args:tt)*],$($t:tt)*)=>{$crate::define_enum_scan!(@field_ty@$tag[$($args)*]$($t)*)};(@field_ty@$tag:ident[[$($args:tt)*]])=>{$($args)*};(@field_ty@$tag:ident[$($args:tt)*])=>{::std::compile_error!(::std::stringify!($($args)*))};(@field_ty$($t:tt)*)=>{$crate::define_enum_scan!(@field_ty@inner[]$($t)*)};(@tag_expr raw,$iter:ident)=>{ScanSource::next_token($iter)?};(@tag_expr$d:ty,$iter:ident)=>{<$d as Scan>::scan($iter)?};(@variant([$($attr:tt)*]$vis:vis$T:ident$d:tt)[$($vars:tt)*])=>{$crate::define_enum_scan!{@def$($attr)*$vis enum$T:$d{$($vars)*}}};(@variant$ctx:tt[$($vars:tt)*]$p:pat=>$v:ident{$($fs:tt)*}$($rest:tt)*)=>{$crate::define_enum_scan!{@field$ctx[$($vars)*]$p=>$v[]$($fs)*;$($rest)*}};(@variant$ctx:tt[$($vars:tt)*]$p:pat=>$v:ident$($rest:tt)*)=>{$crate::define_enum_scan!{@variant$ctx[$($vars)*$p=>$v,]$($rest)*}};(@variant$ctx:tt[$($vars:tt)*],$($rest:tt)*)=>{$crate::define_enum_scan!{@variant$ctx[$($vars)*]$($rest)*}};(@endfield$ctx:tt[$($vars:tt)*]$p:pat=>$v:ident[$($fs:tt)*][$f:ident:$($spec:tt)*],$($rest:tt)*)=>{$crate::define_enum_scan!{@field$ctx[$($vars)*]$p=>$v[$($fs)*[$f:$($spec)*]]$($rest)*}};(@endfield$ctx:tt[$($vars:tt)*]$p:pat=>$v:ident[$($fs:tt)*][$f:ident:$($spec:tt)*];$($rest:tt)*)=>{$crate::define_enum_scan!{@variant$ctx[$($vars)*$p=>$v{$($fs)*[$f:$($spec)*]},]$($rest)*}};(@field$ctx:tt[$($vars:tt)*]$p:pat=>$v:ident[$($fs:tt)*];$($rest:tt)*)=>{$crate::define_enum_scan!{@variant$ctx[$($vars)*$p=>$v{$($fs)*},]$($rest)*}};(@field$ctx:tt[$($vars:tt)*]$p:pat=>$v:ident[$($fs:tt)*]$f:ident:($($tuple:tt)*)$sep:tt$($rest:tt)*)=>{$crate::define_enum_scan!{@endfield$ctx[$($vars)*]$p=>$v[$($fs)*][$f:($($tuple)*)]$sep$($rest)*}};(@field$ctx:tt[$($vars:tt)*]$p:pat=>$v:ident[$($fs:tt)*]$f:ident:[$($x:tt)*]$sep:tt$($rest:tt)*)=>{$crate::define_enum_scan!{@endfield$ctx[$($vars)*]$p=>$v[$($fs)*][$f:[$($x)*]]$sep$($rest)*}};(@field$ctx:tt[$($vars:tt)*]$p:pat=>$v:ident[$($fs:tt)*]$f:ident:$ty:ty=$e:expr,$($rest:tt)*)=>{$crate::define_enum_scan!{@endfield$ctx[$($vars)*]$p=>$v[$($fs)*][$f:$ty=$e],$($rest)*}};(@field$ctx:tt[$($vars:tt)*]$p:pat=>$v:ident[$($fs:tt)*]$f:ident:$ty:ty;$($rest:tt)*)=>{$crate::define_enum_scan!{@endfield$ctx[$($vars)*]$p=>$v[$($fs)*][$f:$ty];$($rest)*}};(@field$ctx:tt[$($vars:tt)*]$p:pat=>$v:ident[$($fs:tt)*]$f:ident:$ty:ty=$e:expr;$($rest:tt)*)=>{$crate::define_enum_scan!{@endfield$ctx[$($vars)*]$p=>$v[$($fs)*][$f:$ty=$e];$($rest)*}};(@field$ctx:tt[$($vars:tt)*]$p:pat=>$v:ident[$($fs:tt)*]$f:ident:$ty:ty,$($rest:tt)*)=>{$crate::define_enum_scan!{@endfield$ctx[$($vars)*]$p=>$v[$($fs)*][$f:$ty],$($rest)*}};(@def$(#[$attr:meta])*$vis:vis enum$T:ident:$d:tt{$($p:pat=>$v:ident$({$([$f:ident:$($spec:tt)*])*})?,)*})=>{$(#[$attr])*$vis enum$T{$($v$({$($f:$crate::define_enum_scan!(@field_ty$($spec)*)),*})?),*}impl Scan for$T{type Output=Self;fn scan(iter:&mut I)->Option{let tag=$crate::define_enum_scan!(@tag_expr$d,iter);match tag{$($p=>{$($(let$f=$crate::scan_value!((*iter),$($spec)*);)*)?Some($T::$v$({$($f),*})?)}),*_=>None,}}}};($(#[$attr:meta])*$vis:vis enum$T:ident:raw{$($body:tt)*})=>{$crate::define_enum_scan!{@variant([$(#[$attr])*]$vis$T raw)[]$($body)*}};($(#[$attr:meta])*$vis:vis enum$T:ident:$d:ty{$($body:tt)*})=>{$crate::define_enum_scan!{@variant([$(#[$attr])*]$vis$T$d)[]$($body)*}};}#[derive(Debug,Copy,Clone)]pub enum Usize1{}impl Scan for Usize1{type Output=usize;fn scan(iter:&mut I)->Option{::scan(iter)?.checked_sub(1)}}#[derive(Debug,Copy,Clone)]pub struct CharWithBase(pub char);impl MarkedScan for CharWithBase{type Output=usize;fn mscan(self,iter:&mut I)->Option{Some((::scan(iter)?as u8-self.0 as u8)as usize)}}#[derive(Debug,Copy,Clone)]pub enum Chars{}impl Scan for Chars{type Output=Vec;fn scan(iter:&mut I)->Option{Some(iter.next_token()?.chars().collect())}}#[derive(Debug,Copy,Clone)]pub struct CharsWithBase(pub char);impl MarkedScan for CharsWithBase{type Output=Vec;fn mscan(self,iter:&mut I)->Option{Some(iter.next_token()?.chars().map(|c|(c as u8-self.0 as u8)as usize).collect())}}#[derive(Debug,Copy,Clone)]pub enum Byte1{}impl Scan for Byte1{type Output=u8;fn scan(iter:&mut I)->Option{let bytes=iter.next_token()?.as_bytes();assert_eq!(bytes.len(),1);Some(bytes[0])}}#[derive(Debug,Copy,Clone)]pub struct ByteWithBase(pub u8);impl MarkedScan for ByteWithBase{type Output=usize;fn mscan(self,iter:&mut I)->Option{Some((::scan(iter)?as u8-self.0)as usize)}}#[derive(Debug,Copy,Clone)]pub enum Bytes{}impl Scan for Bytes{type Output=Vec;fn scan(iter:&mut I)->Option{Some(iter.next_token()?.bytes().collect())}}#[derive(Debug,Copy,Clone)]pub struct BytesWithBase(pub u8);impl MarkedScan for BytesWithBase{type Output=Vec;fn mscan(self,iter:&mut I)->Option{Some(iter.next_token()?.bytes().map(|c|(c-self.0)as usize).collect())}}#[derive(Debug,Copy,Clone)]pub struct Collect::Output>>where T:Scan,B:FromIterator<::Output>{size:usize,_marker:PhantomData(T,B)>}implCollectwhere T:Scan,B:FromIterator<::Output>{pub fn new(size:usize)->Self{Self{size,_marker:PhantomData}}}implMarkedScan for Collectwhere T:Scan,B:FromIterator<::Output>{type Output=B;fn mscan(self,iter:&mut I)->Option{if self.size==0{iter.skip_whitespace();}repeat_with(| |::scan(iter)).take(self.size).collect()}}#[derive(Debug,Copy,Clone)]pub struct SizedCollect::Output>>where T:Scan,B:FromIterator<::Output>{_marker:PhantomData(T,B)>}implScan for SizedCollectwhere T:Scan,B:FromIterator<::Output>{type Output=B;fn scan(iter:&mut I)->Option{let size=usize::scan(iter)?;if size==0{iter.skip_whitespace();}repeat_with(| |::scan(iter)).take(size).collect()}}#[derive(Debug,Copy,Clone)]pub struct Splittedwhere T:Scan{pat:P,_marker:PhantomDataT>}implSplittedwhere T:Scan{pub fn new(pat:P)->Self{Self{pat,_marker:PhantomData}}}implMarkedScan for Splittedwhere T:Scan{type Output=Vec<::Output>;fn mscan(self,iter:&mut I)->Option{let mut iter=Scanner::new_from_iter(iter.next_token()?.split(self.pat));Some(from_fn(| |::scan(&mut iter)).collect())}}implMarkedScan for Splittedwhere T:Scan{type Output=Vec<::Output>;fn mscan(self,iter:&mut I)->Option{let mut iter=Scanner::new_from_iter(iter.next_token()?.split(self.pat));Some(from_fn(| |::scan(&mut iter)).collect())}}implMarkedScan for F where F:Fn(&str)->Option{type Output=T;fn mscan(self,iter:&mut I)->Option{self(iter.next_token()?)}}} // 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())});}