結果
問題 | No.2166 Paint and Fill |
ユーザー | tko919 |
提出日時 | 2022-12-23 22:02:58 |
言語 | C++17 (gcc 12.3.0 + boost 1.83.0) |
結果 |
RE
|
実行時間 | - |
コード長 | 26,269 bytes |
コンパイル時間 | 5,261 ms |
コンパイル使用メモリ | 280,888 KB |
実行使用メモリ | 136,828 KB |
最終ジャッジ日時 | 2024-11-18 04:35:29 |
合計ジャッジ時間 | 79,230 ms |
ジャッジサーバーID (参考情報) |
judge2 / judge3 |
(要ログイン)
テストケース
テストケース表示入力 | 結果 | 実行時間 実行使用メモリ |
---|---|---|
testcase_00 | AC | 832 ms
104,556 KB |
testcase_01 | AC | 786 ms
31,476 KB |
testcase_02 | AC | 1,686 ms
116,304 KB |
testcase_03 | RE | - |
testcase_04 | RE | - |
testcase_05 | RE | - |
testcase_06 | RE | - |
testcase_07 | RE | - |
testcase_08 | RE | - |
testcase_09 | RE | - |
testcase_10 | RE | - |
testcase_11 | RE | - |
testcase_12 | RE | - |
testcase_13 | RE | - |
testcase_14 | RE | - |
testcase_15 | RE | - |
testcase_16 | RE | - |
testcase_17 | RE | - |
testcase_18 | RE | - |
testcase_19 | RE | - |
testcase_20 | RE | - |
testcase_21 | RE | - |
testcase_22 | RE | - |
testcase_23 | RE | - |
testcase_24 | RE | - |
testcase_25 | AC | 2 ms
6,820 KB |
testcase_26 | AC | 2 ms
6,816 KB |
testcase_27 | AC | 2,413 ms
31,724 KB |
testcase_28 | AC | 3,023 ms
31,476 KB |
testcase_29 | AC | 2,473 ms
31,656 KB |
testcase_30 | AC | 3,679 ms
31,724 KB |
testcase_31 | AC | 3,723 ms
31,852 KB |
testcase_32 | AC | 3,711 ms
31,784 KB |
testcase_33 | AC | 3,715 ms
31,724 KB |
testcase_34 | AC | 3,708 ms
31,848 KB |
testcase_35 | AC | 3,776 ms
31,720 KB |
testcase_36 | AC | 3,755 ms
31,724 KB |
testcase_37 | AC | 3,737 ms
31,852 KB |
testcase_38 | AC | 3,737 ms
31,724 KB |
testcase_39 | AC | 3,765 ms
31,724 KB |
ソースコード
#line 1 "library/Template/template.hpp" #include <bits/stdc++.h> using namespace std; #define rep(i,a,b) for(int i=(int)(a);i<(int)(b);i++) #define ALL(v) (v).begin(),(v).end() using ll=long long int; const int inf = 0x3fffffff; const ll INF = 0x1fffffffffffffff; template<typename T>inline bool chmax(T& a,T b){if(a<b){a=b;return 1;}return 0;} template<typename T>inline bool chmin(T& a,T b){if(a>b){a=b;return 1;}return 0;} #line 2 "library/Utility/fastio.hpp" #include <unistd.h> class FastIO{ static constexpr int L=1<<16; char rdbuf[L]; int rdLeft=0,rdRight=0; inline void reload(){ int len=rdRight-rdLeft; memmove(rdbuf,rdbuf+rdLeft,len); rdLeft=0,rdRight=len; rdRight+=fread(rdbuf+len,1,L-len,stdin); } inline bool skip(){ for(;;){ while(rdLeft!=rdRight and rdbuf[rdLeft]<=' ')rdLeft++; if(rdLeft==rdRight){ reload(); if(rdLeft==rdRight)return false; } else break; } return true; } template<typename T,enable_if_t<is_integral<T>::value,int> =0>inline bool _read(T& x){ if(!skip())return false; if(rdLeft+20>=rdRight)reload(); bool neg=false; if(rdbuf[rdLeft]=='-'){ neg=true; rdLeft++; } x=0; while(rdbuf[rdLeft]>='0' and rdLeft<rdRight){ x=x*10+(neg?-(rdbuf[rdLeft++]^48):(rdbuf[rdLeft++]^48)); } return true; } template<typename T,enable_if_t<is_floating_point<T>::value,int> =0>inline bool _read(T& x){ if(!skip())return false; if(rdLeft+20>=rdRight)reload(); bool neg=false; if(rdbuf[rdLeft]=='-'){ neg=true; rdLeft++; } x=0; while(rdbuf[rdLeft]>='0' and rdbuf[rdLeft]<='9' and rdLeft<rdRight){ x=x*10+(rdbuf[rdLeft++]^48); } if(rdbuf[rdLeft]!='.')return true; rdLeft++; T base=.1; while(rdbuf[rdLeft]>='0' and rdbuf[rdLeft]<='9' and rdLeft<rdRight){ x+=base*(rdbuf[rdLeft++]^48); base*=.1; } if(neg)x=-x; return true; } inline bool _read(char& x){ if(!skip())return false; if(rdLeft+1>=rdRight)reload(); x=rdbuf[rdLeft++]; return true; } inline bool _read(string& x){ if(!skip())return false; for(;;){ int pos=rdLeft; while(pos<rdRight and rdbuf[pos]>' ')pos++; x.append(rdbuf+rdLeft,pos-rdLeft); if(rdLeft==pos)break; rdLeft=pos; if(rdLeft==rdRight)reload(); else break; } return true; } template<typename T>inline bool _read(vector<T>& v){ for(auto& x:v){ if(!_read(x))return false; } return true; } char wtbuf[L],tmp[50]; int wtRight=0; inline void flush(){ fwrite(wtbuf,1,wtRight,stdout); wtRight=0; } inline void _write(const char& x){ if(wtRight>L-32)flush(); wtbuf[wtRight++]=x; } inline void _write(const string& x){ for(auto& c:x)_write(c); } template<typename T,enable_if_t<is_integral<T>::value,int> =0>inline void _write(T x){ if(wtRight>L-32)flush(); if(x==0){ _write('0'); return; } else if(x<0){ _write('-'); if (__builtin_expect(x == std::numeric_limits<T>::min(), 0)) { switch (sizeof(x)) { case 2: _write("32768"); return; case 4: _write("2147483648"); return; case 8: _write("9223372036854775808"); return; } } x=-x; } int pos=0; while(x!=0){ tmp[pos++]=char((x%10)|48); x/=10; } rep(i,0,pos)wtbuf[wtRight+i]=tmp[pos-1-i]; wtRight+=pos; } template<typename T>inline void _write(const vector<T>& v){ rep(i,0,v.size()){ if(i)_write(' '); _write(v[i]); } } public: FastIO(){} ~FastIO(){flush();} inline void read(){} template <typename Head, typename... Tail>inline void read(Head& head,Tail&... tail){ assert(_read(head)); read(tail...); } template<bool ln=true,bool space=false>inline void write(){if(ln)_write('\n');} template <bool ln=true,bool space=false,typename Head, typename... Tail>inline void write(const Head& head,const Tail&... tail){ if(space)_write(' '); _write(head); write<ln,true>(tail...); } }; /** * @brief Fast IO */ #line 3 "sol.cpp" #line 2 "library/Math/modint.hpp" template<int mod=1000000007>struct fp { int v; static int get_mod(){return mod;} int inv() const{ int tmp,a=v,b=mod,x=1,y=0; while(b)tmp=a/b,a-=tmp*b,swap(a,b),x-=tmp*y,swap(x,y); if(x<0){x+=mod;} return x; } fp(ll x=0){init(x%mod+mod);} fp& init(ll x){v=(x<mod?x:x-mod); return *this;} fp operator-()const{return fp()-*this;} fp pow(ll t){assert(t>=0); fp res=1,b=*this; while(t){if(t&1)res*=b;b*=b;t>>=1;} return res;} fp& operator+=(const fp& x){return init(v+x.v);} fp& operator-=(const fp& x){return init(v+mod-x.v);} fp& operator*=(const fp& x){v=ll(v)*x.v%mod; return *this;} fp& operator/=(const fp& x){v=ll(v)*x.inv()%mod; return *this;} fp operator+(const fp& x)const{return fp(*this)+=x;} fp operator-(const fp& x)const{return fp(*this)-=x;} fp operator*(const fp& x)const{return fp(*this)*=x;} fp operator/(const fp& x)const{return fp(*this)/=x;} bool operator==(const fp& x)const{return v==x.v;} bool operator!=(const fp& x)const{return v!=x.v;} friend istream& operator>>(istream& is,fp& x){return is>>x.v;} friend ostream& operator<<(ostream& os,const fp& x){return os<<x.v;} }; template<typename T>struct factorial { vector<T> Fact,Finv,Inv; factorial(int maxx){ Fact.resize(maxx); Finv.resize(maxx); Inv.resize(maxx); Fact[0]=Fact[1]=Finv[0]=Finv[1]=Inv[1]=1; rep(i,2,maxx){Fact[i]=Fact[i-1]*i;} Finv[maxx-1]=Fact[maxx-1].inv(); for(int i=maxx-1;i>=2;i--){Finv[i-1]=Finv[i]*i; Inv[i]=Finv[i]*Fact[i-1];} } T fact(int n,bool inv=0){if(n<0)return 0; return (inv?Finv[n]:Fact[n]);} T inv(int n){if(n<0)return 0; return Inv[n];} T nPr(int n,int r,bool inv=0){if(n<0||n<r||r<0)return 0; return fact(n,inv)*fact(n-r,inv^1);} T nCr(int n,int r,bool inv=0){if(n<0||n<r||r<0)return 0; return fact(n,inv)*fact(r,inv^1)*fact(n-r,inv^1);} T nHr(int n,int r,bool inv=0){return nCr(n+r-1,r,inv);} }; /** * @brief Modint */ #line 2 "library/Convolution/ntt.hpp" template<typename T,unsigned p=3>struct NTT{ vector<T> rt,irt; NTT(int lg=21){ unsigned m=T::get_mod()-1; T prt=p; rt.resize(lg); irt.resize(lg); rep(k,0,lg){ rt[k]=-prt.pow(m>>(k+2)); irt[k]=rt[k].inv(); } } void ntt(vector<T>& f,bool inv=0){ int n=f.size(); if(inv){ for(int m=1;m<n;m<<=1){ T w=1; for(int s=0,t=0;s<n;s+=m*2){ for(int i=s,j=s+m;i<s+m;i++,j++){ auto x=f[i],y=f[j]; f[i]=x+y; f[j]=(x-y)*w; } w*=irt[__builtin_ctz(++t)]; } } T mul=T(n).inv(); rep(i,0,n)f[i]*=mul; }else{ for(int m=n;m>>=1;){ T w=1; for(int s=0,t=0;s<n;s+=m*2){ for(int i=s,j=s+m;i<s+m;i++,j++){ auto x=f[i],y=f[j]*w; f[i]=x+y; f[j]=x-y; } w*=rt[__builtin_ctz(++t)]; } } } } vector<T> mult(const vector<T>& a,const vector<T>& b,bool same=0){ if(a.empty() or b.empty())return vector<T>(); int n=a.size()+b.size()-1,m=1<<__lg(n*2-1); vector<T> res(m); rep(i,0,a.size()){res[i]=a[i];} ntt(res); if(same)rep(i,0,m)res[i]*=res[i]; else{ vector<T> c(m); rep(i,0,b.size())c[i]=b[i]; ntt(c); rep(i,0,m)res[i]*=c[i]; } ntt(res,1); res.resize(n); return res; } }; /** * @brief Number Theoretic Transform */ #line 2 "library/FPS/fps.hpp" template<typename T>struct Poly:vector<T>{ Poly(int n=0){this->assign(n,T());} Poly(const vector<T>& f){this->assign(ALL(f));} T eval(const T& x){ T res; for(int i=this->size()-1;i>=0;i--)res*=x,res+=this->at(i); return res; } Poly rev()const{Poly res=*this; reverse(ALL(res)); return res;} void shrink(){while(!this->empty() and this->back()==0)this->pop_back();} vector<T> mult(const vector<T>& a,const vector<T>& b,bool same=0)const{ if(a.empty() or b.empty())return vector<T>(); int n=a.size()+b.size()-1,m=1<<__lg(n*2-1); vector<T> res(m); rep(i,0,a.size())res[i]=a[i]; NTT(res,0); if(same)rep(i,0,m)res[i]*=res[i]; else{ vector<T> c(m); rep(i,0,b.size())c[i]=b[i]; NTT(c,0); rep(i,0,m)res[i]*=c[i]; } NTT(res,1); res.resize(n); return res; } Poly square()const{return Poly(mult(*this,*this,1));} Poly operator-()const{return Poly()-*this;} Poly operator+(const Poly& g)const{return Poly(*this)+=g;} Poly operator+(const T& g)const{return Poly(*this)+=g;} Poly operator-(const Poly& g)const{return Poly(*this)-=g;} Poly operator-(const T& g)const{return Poly(*this)-=g;} Poly operator*(const Poly& g)const{return Poly(*this)*=g;} Poly operator*(const T& g)const{return Poly(*this)*=g;} Poly operator/(const Poly& g)const{return Poly(*this)/=g;} Poly operator%(const Poly& g)const{return Poly(*this)%=g;} Poly& operator+=(const Poly& g){ if(g.size()>this->size())this->resize(g.size()); rep(i,0,g.size()){(*this)[i]+=g[i];} return *this; } Poly& operator+=(const T& g){ if(this->empty())this->push_back(0); (*this)[0]+=g; return *this; } Poly& operator-=(const Poly& g){ if(g.size()>this->size())this->resize(g.size()); rep(i,0,g.size()){(*this)[i]-=g[i];} return *this; } Poly& operator-=(const T& g){ if(this->empty())this->push_back(0); (*this)[0]-=g; return *this; } Poly& operator*=(const Poly& g){ *this=mult(*this,g,0); return *this; } Poly& operator*=(const T& g){ rep(i,0,this->size())(*this)[i]*=g; return *this; } Poly& operator/=(const Poly& g){ if(g.size()>this->size()){ this->clear(); return *this; } Poly g2=g; reverse(ALL(*this)); reverse(ALL(g2)); int n=this->size()-g2.size()+1; this->resize(n); g2.resize(n); *this*=g2.inv(); this->resize(n); reverse(ALL(*this)); shrink(); return *this; } Poly& operator%=(const Poly& g){*this-=*this/g*g; shrink(); return *this;} Poly diff()const{ Poly res(this->size()-1); rep(i,0,res.size())res[i]=(*this)[i+1]*(i+1); return res; } Poly inte()const{ Poly res(this->size()+1); for(int i=res.size()-1;i;i--)res[i]=(*this)[i-1]/i; return res; } Poly log()const{ assert(this->front()==1); const int n=this->size(); Poly res=diff()*inv(); res=res.inte(); res.resize(n); return res; } Poly shift(const int& c)const{ const int n=this->size(); Poly res=*this,g(n); g[0]=1; rep(i,1,n)g[i]=g[i-1]*c/i; vector<T> fact(n,1); rep(i,0,n){ if(i)fact[i]=fact[i-1]*i; res[i]*=fact[i]; } res=res.rev(); res*=g; res.resize(n); res=res.rev(); rep(i,0,n)res[i]/=fact[i]; return res; } Poly inv()const{ const int n=this->size(); Poly res(1); res.front()=T(1)/this->front(); for(int k=1;k<n;k<<=1){ Poly f(k*2),g(k*2); rep(i,0,min(n,k*2))f[i]=(*this)[i]; rep(i,0,k)g[i]=res[i]; NTT(f,0); NTT(g,0); rep(i,0,k*2)f[i]*=g[i]; NTT(f,1); rep(i,0,k){f[i]=0; f[i+k]=-f[i+k];} NTT(f,0); rep(i,0,k*2)f[i]*=g[i]; NTT(f,1); rep(i,0,k)f[i]=res[i]; swap(res,f); } res.resize(n); return res; } Poly exp()const{ const int n=this->size(); if(n==1)return Poly({T(1)}); Poly b(2),c(1),z1,z2(2); b[0]=c[0]=z2[0]=z2[1]=1; b[1]=(*this)[1]; for(int k=2;k<n;k<<=1){ Poly y=b; y.resize(k*2); NTT(y,0); z1=z2; Poly z(k); rep(i,0,k)z[i]=y[i]*z1[i]; NTT(z,1); rep(i,0,k>>1)z[i]=0; NTT(z,0); rep(i,0,k)z[i]*=-z1[i]; NTT(z,1); c.insert(c.end(),z.begin()+(k>>1),z.end()); z2=c; z2.resize(k*2); NTT(z2,0); Poly x=*this; x.resize(k); x=x.diff();x.resize(k); NTT(x,0); rep(i,0,k)x[i]*=y[i]; NTT(x,1); Poly bb=b.diff(); rep(i,0,k-1)x[i]-=bb[i]; x.resize(k*2); rep(i,0,k-1){x[k+i]=x[i]; x[i]=0;} NTT(x,0); rep(i,0,k*2)x[i]*=z2[i]; NTT(x,1); x.pop_back(); x=x.inte(); rep(i,k,min(n,k*2))x[i]+=(*this)[i]; rep(i,0,k)x[i]=0; NTT(x,0); rep(i,0,k*2)x[i]*=y[i]; NTT(x,1); b.insert(b.end(),x.begin()+k,x.end()); } b.resize(n); return b; } Poly pow(ll t){ if(t==0){ Poly res(this->size()); res[0]=1; return res; } int n=this->size(),k=0; while(k<n and (*this)[k]==0)k++; Poly res(n); if(__int128_t(t)*k>=n)return res; n-=t*k; Poly g(n); T c=(*this)[k],ic=c.inv(); rep(i,0,n)g[i]=(*this)[i+k]*ic; g=g.log(); for(auto& x:g)x*=t; g=g.exp(); c=c.pow(t); rep(i,0,n)res[i+t*k]=g[i]*c; return res; } void NTT(vector<T>& a,bool inv)const; }; /** * @brief Formal Power Series (NTT-friendly mod) */ #line 7 "sol.cpp" using Fp=fp<998244353>; NTT<Fp,3> ntt; template<>void Poly<Fp>::NTT(vector<Fp>& v,bool inv)const{return ntt.ntt(v,inv);} #line 2 "library/FPS/samplepointshift.hpp" template<typename T>Poly<T> SamplePointsShift(vector<T>& ys,T c,int m=-1){ ll n=ys.size()-1,C=c.v%T::get_mod(); if(m==-1)m=n+1; factorial<T> fact(ys.size()); if(C<=n){ Poly<T> res; rep(i,C,n+1)res.push_back(ys[i]); if(int(res.size())>=m){ res.resize(m); return res; } auto add=SamplePointsShift<T>(ys,n+1,m-res.size()); for(int i=0;int(res.size())<m;i++){ res.push_back(add[i]); } return res; } if(C+m>T::get_mod()){ auto res=SamplePointsShift<T>(ys,c,T::get_mod()-c.v); auto add=SamplePointsShift<T>(ys,0,m-res.size()); rep(i,0,add.size())res.push_back(add[i]); return res; } Poly<T> A(n+1),B(m+n); rep(i,0,n+1){ A[i]=ys[i]*fact.fact(i,1)*fact.fact(n-i,1); if((n-i)&1)A[i]=-A[i]; } rep(i,0,m+n)B[i]=Fp(1)/(c-n+i); auto AB=A*B; vector<Fp> res(m); Fp base=1; rep(x,0,n+1)base*=(c-x); rep(i,0,m){ res[i]=AB[n+i]*base; base*=(c+i+1); base*=B[i]; } return res; } /** * @brief Shift of Sampling Points of Polynomial */ #line 2 "library/Math/matrix.hpp" template<class T>struct Matrix{ int h,w; vector<vector<T>> val; T det; Matrix(){} Matrix(int n):h(n),w(n),val(vector<vector<T>>(n,vector<T>(n))){} Matrix(int n,int m):h(n),w(m),val(vector<vector<T>>(n,vector<T>(m))){} vector<T>& operator[](const int i){return val[i];} Matrix& operator+=(const Matrix& m){ assert(h==m.h and w==m.w); rep(i,0,h)rep(j,0,w)val[i][j]+=m.val[i][j]; return *this; } Matrix& operator-=(const Matrix& m){ assert(h==m.h and w==m.w); rep(i,0,h)rep(j,0,w)val[i][j]-=m.val[i][j]; return *this; } Matrix& operator*=(const Matrix& m){ assert(w==m.h); Matrix<T> res(h,m.w); rep(i,0,h)rep(j,0,m.w)rep(k,0,w)res.val[i][j]+=val[i][k]*m.val[k][j]; *this=res; return *this; } Matrix operator+(const Matrix& m)const{return Matrix(*this)+=m;} Matrix operator-(const Matrix& m)const{return Matrix(*this)-=m;} Matrix operator*(const Matrix& m)const{return Matrix(*this)*=m;} Matrix pow(ll k){ Matrix<T> res(h,h),c=*this; rep(i,0,h)res.val[i][i]=1; while(k){if(k&1)res*=c; c*=c; k>>=1;} return res; } vector<int> gauss(int c=-1){ if(val.empty())return {}; if(c==-1)c=w; int cur=0; vector<int> res; det=1; rep(i,0,c){ if(cur==h)break; rep(j,cur,h)if(val[j][i]!=0){ swap(val[cur],val[j]); if(cur!=j)det*=-1; break; } det*=val[cur][i]; if(val[cur][i]==0)continue; rep(j,0,h)if(j!=cur){ T z=val[j][i]/val[cur][i]; rep(k,i,w)val[j][k]-=val[cur][k]*z; } res.push_back(i); cur++; } return res; } Matrix inv(){ assert(h==w); Matrix base(h,h*2),res(h,h); rep(i,0,h)rep(j,0,h)base[i][j]=val[i][j]; rep(i,0,h)base[i][h+i]=1; base.gauss(h); rep(i,0,h)rep(j,0,h)res[i][j]=base[i][h+j]/base[i][i]; return res; } bool operator==(const Matrix& m){ assert(h==m.h and w==m.w); rep(i,0,h)rep(j,0,w)if(val[i][j]!=m.val[i][j])return false; return true; } bool operator!=(const Matrix& m){ assert(h==m.h and w==m.w); rep(i,0,h)rep(j,0,w)if(val[i][j]==m.val[i][j])return false; return true; } friend istream& operator>>(istream& is,Matrix& m){ rep(i,0,m.h)rep(j,0,m.w)is>>m[i][j]; return is; } friend ostream& operator<<(ostream& os,Matrix& m){ rep(i,0,m.h){ rep(j,0,m.w)os<<m[i][j]<<(j==m.w-1 and i!=m.h-1?'\n':' '); } return os; } }; /** * @brief Matrix */ #line 3 "library/Math/linearequation.hpp" template<typename T>pair<vector<T>,Matrix<T>> LinearEquation(Matrix<T> a,vector<T> b){ int h=a.h,w=a.w; rep(i,0,h)a[i].push_back(b[i]); a.w++; vector<int> idx=a.gauss(w); rep(i,idx.size(),h)if(a[i][w]!=0)return {{},{}}; vector<T> res(w); rep(i,0,idx.size())res[idx[i]]=a[i][w]/a[i][idx[i]]; Matrix<T> d(w,h+w); rep(i,0,h)rep(j,0,w)d[j][i]=a[i][j]; rep(i,0,w)d[i][h+i]=1; int r=d.gauss(h).size(); Matrix<T> basis(w-r,w); rep(i,r,w)basis[i-r]={d[i].begin()+h,d[i].end()}; return {res,basis}; } /** * @brief Linear Equation */ #line 5 "library/FPS/p-recursive.hpp" template<typename T>Matrix<T> PrefixProdOfPolyMatrix(Matrix<Poly<T>>& m,ll K){ using Mat=Matrix<T>; int n=m.val.size(); int deg=1; rep(i,0,n)rep(j,0,n)chmax(deg,(int)m[i][j].size()-1); ll SQ=1; while(SQ*SQ*deg<K)SQ<<=1; T iSQ=T(SQ).inv(); vector<Mat> G(deg+1); rep(k,0,deg+1){ G[k]=Mat(n,n); rep(i,0,n)rep(j,0,n)G[k][i][j]=m[i][j].eval(SQ*k); } auto process=[&](vector<Mat>& base,T x)->vector<Mat>{ int D=base.size(); vector ret(D,Mat(n,n)); rep(i,0,n)rep(j,0,n){ vector<T> val(D); rep(k,0,D)val[k]=base[k][i][j]; auto add=SamplePointsShift<T>(val,x); rep(k,0,D)ret[k][i][j]=add[k]; } return ret; }; for(ll w=1;w<SQ;w<<=1){ auto G1=process(G,iSQ*w); auto G2=process(G,w*deg+1); auto G3=process(G,iSQ*w+w*deg+1); rep(i,0,w*deg+1)G1[i]*=G[i],G3[i]*=G2[i]; G1.insert(G1.end(),ALL(G3)); G1.pop_back(); swap(G,G1); } Mat ret(n,n); rep(i,0,n)ret[i][i]=1; ll k=0; while(k*SQ+SQ<=K)ret=G[k++]*ret; k*=SQ; while(k<K){ Mat mul(n,n); rep(i,0,n)rep(j,0,n)mul[i][j]=m[i][j].eval(k); ret=mul*ret; k++; } return ret; } // a_{n+i}*f_n(i)+...+a_i*f_0(i)=0 // {f_r}:dec order!!! template<typename T>vector<Poly<T>> FindPRecursive(vector<T>& a,int d){ int n=a.size(); int k=(n+2)/(d+2)-1; if(k<=0)return {}; int m=(d+1)*(k+1); Matrix<T> mat(m-1,m); rep(i,0,m-1)rep(j,0,k+1){ T base=1; rep(deg,0,d+1){ mat[i][(d+1)*j+deg]=a[i+j]*base; base*=(i+j); } } auto basis=LinearEquation(mat,vector<T>(m-1)).second; if(basis.val.empty())return {}; auto c=basis[0]; vector<Poly<T>> ret; for(int i=0;i*(d+1)<(int)c.size();i++){ Poly<T> add,base({T(i),T(1)}); for(int j=d;j>=0;j--){ add*=base; if(c[i*(d+1)+j]!=0)add+=c[i*(d+1)+j]; } ret.push_back(add); } while(ret.back().empty())ret.pop_back(); reverse(ALL(ret)); return ret; } template<typename T>T KthtermOfPRecursive(vector<T>& a,vector<Poly<T>>& fs,ll k){ int n=fs.size()-1; assert(int(a.size())>=n); if(k<int(a.size()))return a[k]; Matrix<Poly<T>> m(n),den(1); Matrix<T> base(n); rep(i,0,n)m[0][i]=-fs[i+1]; rep(i,1,n)m[i][i-1]=fs[0]; den[0][0]=fs[0]; rep(i,0,n)base[i][0]=a[n-1-i]; T ret=(PrefixProdOfPolyMatrix(m,k-n+1)*base)[0][0]; ret/=PrefixProdOfPolyMatrix(den,k-n+1)[0][0]; return ret; } template<typename T>T KthtermEsper(vector<T>& a,ll k){ if(k<(int)a.size())return a[k]; int n=a.size()-1; vector<Fp> b=a; b.pop_back(); for(int d=0;;d++){ if((n+2)/(d+2)<=1)break; auto fs=FindPRecursive(b,d); if(KthtermOfPRecursive(b,fs,n)==a.back()){ return KthtermOfPRecursive(a,fs,k); } } cerr<<"esper Failed"<<'\n'; assert(0); } /** * @brief P-recursive */ #line 2 "library/FPS/multieval.hpp" template<typename T>struct MultiEval{ int m,n; vector<Poly<T>> t; MultiEval(vector<T>& v){ m=v.size(),n=1; while(n<m)n<<=1; t.resize(n<<1); rep(i,0,n){ T w=(i<m?v[i]:0); t[n+i]=Poly<T>({-w,T(1)}); } for(int i=n-1;i;i--)t[i]=t[i*2]*t[i*2+1]; } vector<T> run(const vector<T>& f){ vector<Poly<T>> c(n*2); auto v=t[1].rev(); v.resize(f.size()); v=v.inv().rev()*Poly<T>(f); v.erase(v.begin(),v.begin()+f.size()-1); v.resize(n); reverse(ALL(v)); c[1]=v; rep(i,1,n){ int d=c[i].size(); rep(k,0,2){ auto add=t[i*2+(k^1)]; add.resize(d/2+1); add=add.rev(); add*=c[i]; add.resize(d); c[i*2+k]=vector<T>(add.begin()+d/2,add.end()); } } vector<T> res(m); rep(i,0,m)res[i]=c[n+i][0]; return res; } vector<T> build(vector<T>& ys){ auto w=t[1].rev(); w.resize(m+1); auto vs=run(w.rev().diff()); rep(i,0,m)ys[i]/=vs[i]; vector<Poly<T>> c(n*2); rep(i,0,n){ if(i<m)c[n+i]=Poly<T>({ys[i]}); else c[n+i]=Poly<T>({T()}); } for(int i=n-1;i;i--)c[i]=c[i*2]*t[i*2+1]+c[i*2+1]*t[i*2]; c[1]=vector<T>(c[1].begin()+(n-m),c[1].end()); c[1].resize(m); return c[1]; } }; /** * @brief Multipoint Evaluation */ #line 13 "sol.cpp" FastIO io; void solve1(int t){ vector<ll> n(t),k(t); rep(i,0,t)io.read(n[i],k[i]); int m=1<<17; vector<Poly<Fp>> subprod(m*2,Poly<Fp>({Fp(1)})); using P=pair<ll,ll>; vector que(m,vector<P>()); rep(i,0,t){ que[k[i]].push_back({n[i],i}); } rep(k,0,m)if(que[k].size()){ deque<Poly<Fp>> deq; for(auto& [N,_]:que[k])deq.push_back(Poly<Fp>({Fp(-N),Fp(1)})); while(deq.size()>1){ auto A=deq.front(); deq.pop_front(); auto B=deq.front(); deq.pop_front(); deq.push_back(A*B); } subprod[m+k]=deq.front(); } for(int i=m-1;i;i--)subprod[i]=subprod[i*2]*subprod[i*2+1]; vector<Fp> ret(t); vector mat(m*2,Matrix<Poly<Fp>>(2)); auto dfs=[&](auto& dfs,int L,int R,int id)->void{ if(R-L==1){ if(que[L].size()){ vector<Fp> xs; for(auto& [x,_]:que[L])xs.push_back(x); MultiEval<Fp> buf(xs); auto ys=buf.run(mat[id][0][0]); rep(i,0,que[L].size())ret[que[L][i].second]=ys[i]; } mat[id][0][0]=Poly<Fp>({Fp(-2*L),Fp(2)}); mat[id][0][1]=Poly<Fp>({Fp(-L)*(L-1)/2,Fp(L)}); mat[id][1][0]=Poly<Fp>({Fp(1)}); mat[id][1][1]={}; return; } int mid=(L+R)>>1; mat[id*2]=mat[id]; rep(i,0,2)rep(j,0,2)mat[id*2][i][j]%=subprod[id*2]; dfs(dfs,L,mid,id*2); mat[id*2+1]=mat[id*2]*mat[id]; rep(i,0,2)rep(j,0,2)mat[id*2+1][i][j]%=subprod[id*2+1]; dfs(dfs,mid,R,id*2+1); mat[id]=mat[id*2+1]*mat[id*2]; rep(i,0,2)rep(j,0,2)mat[id][i][j]%=subprod[id]; return; }; mat[1][0][0]=mat[1][1][1]=Poly<Fp>({Fp(1)}); dfs(dfs,0,m,1); rep(i,0,t)io.write(ret[i].v); } void solve2(int t){ while(t--){ ll n,k; io.read(n,k); if(k>=Fp::get_mod())io.write(0); else{ vector<Fp> a(2); a[0]=1,a[1]=n*2; vector<Poly<Fp>> fs(3); fs[0]=Poly<Fp>({Fp(1)}); fs[1]=Poly<Fp>({-n*2+2,2}); fs[2]=Poly<Fp>({-n,Fp(1-n*2)/2,Fp(1)/2}); Fp ret=KthtermOfPRecursive(a,fs,k); io.write(ret.v); } } } int main(){ int t; io.read(t); if(t>5)solve1(t); else solve2(t); return 0; }