#include #include using namespace std; using mint=atcoder::modint998244353; template struct Number_Theoretic_Transform { static vectordw,dw_inv; static int log; static mint root; static void ntt(vector& f) { init(); const int n=f.size(); for(int m=n;m>>=1;) { mint w=1; for(int s=0,k=0;s& f, bool flag=true) { init(); const int n=f.size(); for(int m=1;m>=1; log++; } dw.resize(log); dw_inv.resize(log); for(int i=0;i>(i+2)); dw_inv[i]=dw[i].inv(); } } }; template vectorNumber_Theoretic_Transform::dw=vector(); template vectorNumber_Theoretic_Transform::dw_inv=vector(); template int Number_Theoretic_Transform::log=0; template mint Number_Theoretic_Transform::root=mint(3); template struct Formal_Power_Series:vector { using FPS=Formal_Power_Series; using vector::vector; using NTT=Number_Theoretic_Transform; void ntt(){NTT::ntt(*this);} void intt(bool flag=true){NTT::intt(*this,flag);} Formal_Power_Series(const vector&x, const vector&y) { const int n=x.size(); auto rec=[&](const auto&rec, int l, int r)->FPS { if(l+1==r)return FPS{-x[l],1}; int m=(l+r)/2; return rec(rec,l,m)*rec(rec,m,r); }; FPS g=rec(rec,0,n); vector dg=g.diff().multipoint_evaluation(x); auto rec2=[&](const auto&rec2, int l, int r)->pair { if(l+1==r)return make_pair(FPS{y[l]/dg[l]},FPS{-x[l],1}); int m=(l+r)/2; auto[fl,gl]=rec2(rec2,l,m); auto[fr,gr]=rec2(rec2,m,r); return make_pair(fl*gr+fr*gl,gl*gr); }; *this=rec2(rec2,0,n).first; } FPS &operator+=(const mint& r) { if(this->empty())this->resize(1); (*this)[0]+=r; return *this; } FPS &operator-=(const mint& r) { if(this->empty())this->resize(1); (*this)[0]-=r; return *this; } FPS &operator*=(const mint& r) { for(mint &x:*this)x*=r; return *this; } FPS &operator/=(const mint& r) { mint invr=r.inv(); for(mint &x:*this)x*=invr; return *this; } FPS operator+(const mint& r)const{return FPS(*this)+=r;} FPS operator-(const mint& r)const{return FPS(*this)-=r;} FPS operator*(const mint& r)const{return FPS(*this)*=r;} FPS operator/(const mint& r)const{return FPS(*this)/=r;} FPS& operator+=(const FPS& f) { if(this->size()resize(f.size()); for(int i=0;i<(int)f.size();i++)(*this)[i]+=f[i]; return *this; } FPS& operator-=(const FPS& f) { if(this->size()resize(f.size()); for(int i=0;i<(int)f.size();i++)(*this)[i]-=f[i]; return *this; } FPS& operator*=(const FPS& f) { *this=convolution(*this,f); return *this; } FPS& operator/=(const FPS& f) { return *this*=f.inv(); } FPS& operator%=(const FPS& f) { *this-=this->div(f)*f; this->shrink(); return *this; } FPS operator+(const FPS& f)const{return FPS(*this)+=f;} FPS operator-(const FPS& f)const{return FPS(*this)-=f;} FPS operator*(const FPS& f)const{return FPS(*this)*=f;} FPS operator/(const FPS& f)const{return FPS(*this)/=f;} FPS operator%(const FPS& f)const{return FPS(*this)%=f;} FPS operator-()const { FPS res(this->size()); for(int i=0;i<(int)this->size();i++)res[i]-=(*this)[i]; return res; } FPS div(FPS f)const { if(this->size()size()-f.size()+1; return (rev().pre(n)*f.rev().inv(n)).pre(n).rev(n); } FPS pre(int deg)const { return FPS(begin(*this),begin(*this)+min((int)this->size(),deg)); } FPS rev(int deg=-1)const { FPS res(*this); if(deg!=-1)res.resize(deg,0); reverse(begin(res),end(res)); return res; } void shrink() { while(!this->empty()&&this->back()==0)this->pop_back(); } FPS dot(FPS f)const { int n=min(this->size(),f.size()); FPS res(n); for(int i=0;i>(int deg)const { if((int)this->size()<=deg)return{}; FPS res(*this); res.erase(res.begin(),res.begin()+deg); return res; } FPS& operator>>=(int deg) { return *this=*this>>(deg); } mint operator()(const mint& r)const { mint res=0,powr=1; for(auto x:*this) { res+=x*powr; powr*=r; } return res; } FPS diff()const { int n=this->size(); FPS res(max(0,n-1)); for(int i=1;isize(); FPS res(n+1); res[0]=0; for(int i=0;isize(); if(deg==-1)deg=n; FPS res(deg); res[0]={(*this)[0].inv()}; for(int d=1;dsize(); vectorinv; inv.reserve(deg+1); inv.push_back(mint::raw(0)); inv.push_back(mint::raw(1)); auto inplace_integral=[&](FPS& f)->void { int n=f.size(); long long mod=mint::mod(); while(inv.size()<=f.size()) { int i=inv.size(); inv.push_back((-inv[mod%i])*(mod/i)); } f.insert(begin(f),mint::raw(0)); for(int i=1;i<=n;i++)f[i]*=inv[i]; }; auto inplace_diff=[](FPS& f)->void { if(f.empty())return; f.erase(begin(f)); mint cef=1; for(int i=0;i<(int)f.size();i++) { f[i]*=cef; cef++; } }; FPS b={1,1size()?(*this)[1]:0}; FPS c={1},z1,z2={1,1}; for(int m=2;msize()),m)); inplace_diff(x); x.push_back(mint::raw(0)); x.ntt(); x=x.dot(y); x.intt(); x-=b.diff(); x.resize(2*m); for(int i=0;isize()),2*m);i++)x[i]+=(*this)[i]; fill(begin(x),begin(x)+m,mint::raw(0)); x.ntt(); x=x.dot(y); x.intt(); b.insert(end(b),begin(x)+m,end(x)); } return FPS{begin(b),begin(b)+deg}; } FPS log(int deg=-1)const { assert((*this)[0]==1); int n=this->size(); if(deg==-1)deg=n; return (this->diff()*this->inv()).pre(deg-1).integral(); } FPS pow(long long k, int deg=-1)const { if(deg==-1)deg=this->size(); if(k==0) { FPS res(deg); res[0]=mint::raw(1); return res; } FPS res=*this; int cnt0=0; while(cnt0(deg-1)/k) { FPS res(deg); return res; } res=res>>cnt0; deg-=cnt0*k; res=((res/res[0]).log(deg)*k).exp(deg)*res[0].pow(k); res=res<<(cnt0*k); return res.pre(deg); } FPS sqrt(int deg=-1)const { auto sqrt_mod=[](mint r)->mint { const int mod=mint::mod(); if(r==0||r==1)return r; if(r.pow((mod-1)>>1)!=1)return -1; mint b=1; while(b.pow((mod-1)>>1)==1)b++; int m=mod-1,e=0; while((m&1)==0)m>>=1,e++; mint x=r.pow((m-1)>>1); mint y=r*x*x; x*=r; mint z=b.pow(m); while(y!=1) { int j=0; mint t=y; while(t!=1)j++,t*=t; z=z.pow(1LL<<(e-j-1)); x*=z; z*=z; y*=z; e=j; } return x; }; if(deg==-1)deg=this->size(); int cnt0=0; while(cnt0size()&&(*this)[cnt0]==0)cnt0++; if(cnt0) { if(cnt0==(int)this->size())return FPS(deg); if(cnt0&1)return{}; if(2*deg<=cnt0)return FPS(deg); FPS res=(*this>>cnt0).sqrt(deg-cnt0/2); if(res.empty())return{}; res=res<<(cnt0/2); res.resize(deg,0); return res; } mint sqr=sqrt_mod((*this)[0]); if(sqr*sqr!=(*this)[0])return{}; FPS res={sqr}; mint inv2=mint(2).inv(); for(int i=1;ipre(i<<1)*res.inv(i<<1))*inv2; return res.pre(deg); } FPS circular_mod(int m)const { FPS res(m); for(int i=0;i<(int)this->size();i++)res[i%m]+=(*this)[i]; return res; } FPS taylor_shift(mint c)const { int n=this->size(); FPS fact(n),fact_inv(n); { // calc fact and fact inv fact[0]=1; for(int i=1;i=1;i--)fact_inv[i-1]=i*fact_inv[i]; } FPS res(*this); res=res.dot(fact); res=res.rev(); FPS bs(n,mint::raw(1)); for(int i=1;imultipoint_evaluation(const vector&x)const { if(x.empty())return{}; int m=x.size(),n=1; if(this->size()==0){return vector(m,0);} if(this->size()==1){return vector(m,(*this)[0]);} while(m>n)n<<=1; vectorf(n<<1,FPS({mint(1)})); for(int i=0;i0;i--)f[i]=f[i<<1]*f[(i<<1)|1]; f[1]=(*this)%f[1]; for(int i=2;i>1]%f[i]; vectorres(m); for(int i=0;i; vectorBerlekamp_Massey(const vector&a) { int n=a.size(); vectorb,c; b.reserve(n+1);b.push_back(1); c.reserve(n+1);c.push_back(1); mint y=1; for(int d=1;d<=n;d++) { int k=b.size(),l=c.size(); mint x=0; for(int i=0;itmp=c; c.insert(c.begin(),k-l,0); for(int i=0;i=q.size()) { FPS r=p.div(q); p-=q*r; p.shrink(); if(k<(int)r.size())res+=r[k]; } if(p.empty())return res; p.resize(q.size()-1); auto sub=[&](const FPS& f, bool odd=0) { FPS g((f.size()+!odd)/2); for(int i=odd;i<(int)f.size();i+=2)g[i>>1]=f[i]; return g; }; while(k) { FPS q2=q; for(int i=1;i<(int)q2.size();i+=2)q2[i]=-q2[i]; p=sub(p*q2,k&1); q=sub(q*q2); k>>=1; } return res+p[0]; } mint linear_recurrence(FPS a, FPS c, long long k) { assert(a.size()==c.size()); c=FPS{1}-(c<<1); return Bostan_Mori((a*c).pre(a.size()),c,k); } mint BMBM(vector&x, long long k) { vectortmp=Berlekamp_Massey(x); int n=tmp.size()-1; FPS a(n),c(n); for(int i=0;i>N>>K; int M=4*K; vectorF(M); F[0]=F[1]=1; vectorA(M); A[0]=1,A[1]=2; for(int i=2;i