多合一无聊题
mod k=t,并且k是p-1的约数
单位根反演石锤了。
所以直接设f[i]表示走i步的方案数,
然后C(L,i)分配位置,再A^i进行矩乘得到f[i]
变成生成函数F(x)=∑f[i]=(A*x+I)^L
求指数mod k=t的系数的和
偏移之后,进行单位根反演
对于t都要求?
NTT
然后WA了
因为要任意模数NTT,。。。。
然后套用全家桶有了9K的代码:
#include<bits/stdc++.h> #define reg register int #define il inline #define fi first #define se second #define mk(a,b) make_pair(a,b) #define numb (ch^'0') using namespace std; typedef long long ll; template<class T>il void rd(T &x){ char ch;x=0;bool fl=false; while(!isdigit(ch=getchar()))(ch=='-')&&(fl=true); for(x=numb;isdigit(ch=getchar());x=x*10+numb); (fl==true)&&(x=-x); } template<class T>il void output(T x){if(x/10)output(x/10);putchar(x%10+'0');} template<class T>il void ot(T x){if(x<0) putchar('-'),x=-x;output(x);putchar(' ');} template<class T>il void prt(T a[],int st,int nd){for(reg i=st;i<=nd;++i) ot(a[i]);putchar(' ');} //--------------------------------------------------------------------------------------------------------------------// namespace Miracle{ const int N=65536*8; int mod,w; int X,L,Y,K,G,GI; const double Pi=acos(-1); il int qm(int x,ll y){int ret=1;while(y){if(y&1) ret=(ll)ret*x%mod;x=(ll)x*x%mod;y>>=1;}return ret;} il int ad(int x,int y){return x+y>=mod?x+y-mod:x+y;} il int sub(int x,int y){return ad(x,mod-y);} il int mul(int x,int y){return (ll)x*y%mod;} struct tr{ int a[3][3]; tr(){memset(a,0,sizeof a);} void init(){ a[0][0]=a[1][1]=a[2][2]=1; } tr friend operator *(const tr &a,const tr &b){ tr c; for(reg i=0;i<3;++i){ for(reg k=0;k<3;++k){ for(reg j=0;j<3;++j){ c.a[i][j]=ad(c.a[i][j],mul(a.a[i][k],b.a[k][j])); } } }return c; } tr friend operator +(const tr &a,const tr &b){ tr c; for(reg i=0;i<3;++i){ for(reg j=0;j<3;++j){ c.a[i][j]=ad(a.a[i][j],b.a[i][j]); } } return c; } tr friend operator -(const tr &a,const tr &b){ tr c; for(reg i=0;i<3;++i){ for(reg j=0;j<3;++j){ c.a[i][j]=ad(a.a[i][j],mod-b.a[i][j]); } } return c; } tr friend operator *(const tr &a,const int &v){ tr c; for(reg i=0;i<3;++i){ for(reg j=0;j<3;++j){ c.a[i][j]=mul(a.a[i][j],v); } } return c; } void out(){ for(reg i=0;i<3;++i){ for(reg j=0;j<3;++j){ cout<<a[i][j]<<" "; }cout<<endl; } } }S,A,I,B,c[N],ans[N]; tr qm(tr x,int y){ tr ret;ret.init(); while(y){ if(y&1) ret=ret*x; x=x*x; y>>=1; } return ret; } namespace Polynomial{ struct Poly{ vector<int>f; Poly(){f.clear();} il int &operator[](const int &x){return f[x];} il const int &operator[](const int &x) const {return f[x];} il void resize(const int &n){f.resize(n);} il int size() const {return f.size();} il void cpy(Poly &b){f.resize(b.size());for(reg i=0;i<(int)f.size();++i)f[i]=b[i];} il void rev(){reverse(f.begin(),f.end());} il void clear(){f.clear();} il void read(const int &n){f.resize(n);for(reg i=0;i<n;++i)rd(f[i]);} il void out() const {for(reg i=0;i<(int)f.size();++i)ot(f[i]);putchar(' ');} }R; il int init(const int &n){int m;for(m=1;m<n;m<<=1);return m;} template<class T>il void rev(T &f){ int lp=f.size(); if(R.size()!=f.size()) { R.resize(f.size()); for(reg i=0;i<lp;++i){ R[i]=(R[i>>1]>>1)|((i&1)?lp>>1:0); } } for(reg i=0;i<lp;++i){ if(i<R[i]) swap(f[i],f[R[i]]); } } } using namespace Polynomial; //--------------------------------------------------------------------------------------------------------------------// il void operator +=(Poly &f,const Poly &g){for(reg i=0;i<f.size();++i) f[i]=ad(f[i],g[i]);} il void operator +=(Poly &f,const int &c){f[0]=ad(f[0],c);} il Poly operator +(Poly f,const Poly &g){for(reg i=0;i<f.size();++i) f[i]=ad(f[i],g[i]);return f;} il Poly operator +(Poly f,const int &c){f[0]=ad(f[0],c);return f;} il void operator -=(Poly &f,const Poly &g){for(reg i=0;i<f.size();++i) f[i]=sub(f[i],g[i]);} il void operator -=(Poly &f,const int &c){f[0]=sub(f[0],c);} il Poly operator -(Poly f,const Poly &g){for(reg i=0;i<f.size();++i) f[i]=sub(f[i],g[i]);return f;} il Poly operator -(Poly f,const int &c){f[0]=sub(f[0],c);return f;} il Poly operator -(Poly f){for(reg i=0;i<f.size();++i) f[i]=mod-f[i];return f;} //--------------------------------------------------------------------------------------------------------------------// namespace FastFourierTransform{ struct cplx{ double x,y; cplx(){x=0.0;y=0.0;} cplx(double xx,double yy){x=xx;y=yy;} cplx friend operator !(cplx a){return cplx(a.x,-a.y);} cplx friend operator +(cplx a,cplx b){return cplx(a.x+b.x,a.y+b.y);} cplx friend operator -(cplx a,cplx b){return cplx(a.x-b.x,a.y-b.y);} cplx friend operator *(cplx a,cplx b){return cplx(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);} }; struct Cps{ vector<cplx>f; Cps(){f.clear();} il cplx &operator[](const int &x){return f[x];} il const cplx &operator[](const int &x) const {return f[x];} il void resize(const int &n){f.resize(n);} il int size() const {return f.size();} il void cpy(Cps &b){f.resize(b.size());for(reg i=0;i<(int)f.size();++i)f[i]=b[i];} il void rev(){reverse(f.begin(),f.end());} il void clear(){f.clear();} il void out(){ for(reg i=0;i<(int)f.size();++i){ cout<<"("<<f[i].x<<","<<f[i].y<<") "; }cout<<endl; } }W; il void FFT(Cps &f,int c){ int n=f.size();rev(f); for(reg p=2;p<=n;p<<=1){ int len=p/2; for(reg l=0;l<n;l+=p){ for(reg k=l;k<l+len;++k){ cplx tmp=f[k+len]*(c>0?W[n/p*(k-l)]:!W[n/p*(k-l)]); f[k+len]=f[k]-tmp; f[k]=f[k]+tmp; } } } if(c==-1){ for(reg i=0;i<n;++i){ f[i].x/=n;f[i].y/=n; } } } il void prework(int n){ if(W.size()!=n){ W.resize(n); for(reg i=0;i<n;++i){ W[i]=cplx(cos(2*Pi/n*i),sin(2*Pi/n*i)); } } } il Poly MTT(const Poly &F,const Poly &G,const int &P){ int n=F.size(),m=G.size(); Cps a,b,c,d; int len=init(n+m-1); a.resize(len);b.resize(len); c.resize(len);d.resize(len); for(reg i=0;i<n;++i){ a[i].x=F[i]>>15;a[i].y=F[i]&32767; } for(reg i=0;i<m;++i){ b[i].x=G[i]>>15;b[i].y=G[i]&32767; } prework(len); FFT(a,1);FFT(b,1); cplx ka,kb,ba,bb; cplx aaa=cplx(0.5,0),bbb=cplx(0,-0.5),o=cplx(0,1); for(reg i=0;i<len;++i){ int j=(len-i)%len; ka=(a[i]+!a[j])*aaa;ba=(a[i]-!a[j])*bbb; kb=(b[i]+!b[j])*aaa;bb=(b[i]-!b[j])*bbb; c[i]=ka*kb+ba*kb*o; d[i]=bb*ka+bb*ba*o; } FFT(c,-1);FFT(d,-1); Poly ret; ret.resize(n+m-1); for(reg i=0;i<n+m-1;++i){ ll A=(ll)(c[i].x+0.5)%P,B=(ll)(c[i].y+0.5)%P; ll C=(ll)(d[i].x+0.5)%P,D=(ll)(d[i].y+0.5)%P; ret[i]=((((A<<30)%P)+((B+C)<<15)%P)%P+D)%P; } return ret; } il void operator *=(Poly &f,Poly g){ f=MTT(f,g,mod); } il void operator *=(Poly &f,const int &c){for(reg i=0;i<f.size();++i) f[i]=mul(f[i],c);} il Poly operator *(Poly f,const Poly &g){f*=g;return f;} il Poly operator *(Poly f,const int &c){for(reg i=0;i<f.size();++i) f[i]=mul(f[i],c);return f;} il Poly Inv(const Poly &f,int n){ if(n==1){ Poly g;g.resize(1);g[0]=qm(f[0],mod-2);return g; } Poly h=Inv(f,(n+1)>>1); Poly tmp=h,t; t.resize(n); for(reg i=0;i<n;++i) t[i]=f[i]; tmp=tmp*tmp*t; h.resize(tmp.size()); Poly g=h*2-tmp; g.resize(n); return g; } } using namespace FastFourierTransform; int pri[N],cnt; void fin(){ int lp=mod-1; for(reg i=2;(ll)i*i<=lp;++i){ if(lp%i==0){ ++cnt; pri[cnt]=i; while(lp%i==0) lp/=i; } } if(lp) pri[++cnt]=lp; // cout<<" cnt "<<cnt<<endl; // prt(pri,1,cnt); G=2; lp=mod-1; for(;;++G){ bool fl=true; for(reg j=1;j<=cnt;++j){ if(qm(G,lp/pri[j])==1) fl=false; } if(fl) break; } } int main(){ int n; rd(n);rd(K);rd(L);rd(X);rd(Y);rd(mod); --X;--Y; fin(); GI=qm(G,mod-2); // cout<<" GG "<<G<<" GI "<<GI<<endl; for(reg i=0;i<n;++i){ for(reg j=0;j<n;++j){ rd(A.a[i][j]); } } I.init(); S.a[0][X]=1; int now=1; w=qm(G,(mod-1)/K); int T=qm(w,mod-2); for(reg i=0;i<K;++i){ c[i]=qm((A*now)+I,L)*qm(w,(ll)i*(i-1)/2); now=mul(now,w); // cout<<" i "<<i<<endl; // c[i].out(); // cout<<endl; } Poly g; g.resize(2*K-1); for(reg i=0;i<2*K-1;++i){ g[i]=qm(T,(ll)i*(i-1)/2); } g.rev(); Poly f; f.resize(2*K-1); for(reg i=0;i<3;++i){ int j=Y; f.clear(); f.resize(2*K-1); for(reg k=0;k<K;++k){ f[k]=c[k].a[i][j]; } // f.out(); // g.out(); // cout<<endl; f*=g; // cout<<" ff "<<endl; // f.out(); // cout<<" edn "<<endl; for(reg k=0;k<f.size();++k){ ans[k].a[i][j]=f[k]; } // for(reg k=0;k<K;++k){ // // cout<<" ans[k] "<<k<<endl; // // ans[k].out(); // } } for(reg t=0;t<K;++t){ tr now=ans[2*K-2-t]; // cout<<" tt "<<t<<endl;now.out(); now=now*qm(w,(ll)t*(t-1)/2)*qm(K,mod-2); tr out=S*now; // ot(out.a[0][Y]); printf("%d ",out.a[0][Y]); } return 0; } } signed main(){ Miracle::main(); return 0; } /* Author: *Miracle* Date: 2019/4/8 18:57:00 */