[HNOI2019]白兔之舞(矩阵快速幂+单位根反演)
非常抱歉,这篇文章鸽了。真的没时间写了。
#include<bits/stdc++.h> using namespace std; typedef long long ll; #define cp complex<long double> #define pi acosl(-1) const int N=6e5+7; struct mat{ll a[5][5];}A1,A2; ll n,k,l,X,Y,mod,w,nn,ans[N],a1[N],a2[N],R[N]; cp A[N],B[N],C[N],D[N],E[N],F[N],G[N]; ll qpow(ll a,ll b) { ll ret=1; while(b) { if(b&1)ret=ret*a%mod; a=a*a%mod,b>>=1; } return ret; } mat operator*(mat a,mat b) { mat c; for(int i=1;i<=n;i++) for(int j=1;j<=n;j++) { c.a[i][j]=0; for(int k=1;k<=n;k++)c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j])%mod; } return c; } mat qpow(mat a,ll b) { mat ret={0}; for(int i=1;i<=n;i++) for(int j=1;j<=n;j++) if(i==j)ret.a[i][j]=1; while(b) { if(b&1)ret=ret*a; a=a*a,b>>=1; } return ret; } ll getrt(ll x) { vector<int>vec;vec.clear(); for(ll i=2;i*i<=(x-1);i++)if((x-1)%i==0)vec.push_back(i),vec.push_back((x-1)/i); for(ll i=2;;i++) { bool c=1; for(unsigned j=0;j<vec.size()&&c;j++)if(qpow(i,vec[j])==1)c=0; if(c)return i; } } ll c1(ll x){return x*(x-1)/2;} void FFT(cp*a,int f) { for(int i=0;i<nn;i++)if(i<R[i])swap(a[i],a[R[i]]); for(int i=1;i<nn;i*=2) { cp wn(cosl(pi/i),sinl(pi/i)); for(int j=0;j<nn;j+=i*2) { cp w(1,0),x,y; for(int k=0;k<i;k++,w=w*wn) x=a[j+k],y=a[j+k+i]*w,a[j+k]=x+y,a[j+k+i]=x-y; } } if(f==1)return; reverse(a+1,a+nn); for(int i=0;i<nn;i++)a[i]/=nn; } int main() { cin>>n>>k>>l>>X>>Y>>mod; for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)cin>>A1.a[i][j]; w=qpow(getrt(mod),(mod-1)/k); for(int j=0;j<k;j++) { for(int i=1;i<=n;i++) for(int p=1;p<=n;p++) { A2.a[i][p]=A1.a[i][p]*qpow(w,j)%mod; if(i==p)A2.a[i][p]++; } A2=qpow(A2,l),a2[j]=A2.a[X][Y]; } for(int j=0;j<k;j++)a2[j]=a2[j]*qpow(w,c1(j))%mod; for(int j=0;j<=k*2;j++)a1[j]=qpow(w,c1(j)*(mod-2)); reverse(a1,a1+k*2+1); int m=k*3,L=0; for(int i=0;i<m;i++)A[i]=a2[i]/32768,B[i]=a2[i]%32768; for(int i=0;i<m;i++)C[i]=a1[i]/32768,D[i]=a1[i]%32768; for(nn=1;nn<=m;nn<<=1)L++; for(int i=0;i<nn;i++)R[i]=(R[i>>1]>>1)|((i&1)<<L-1); FFT(A,1),FFT(B,1),FFT(C,1),FFT(D,1); for(int i=0;i<nn;i++)E[i]=A[i]*C[i],F[i]=A[i]*D[i]+B[i]*C[i],G[i]=B[i]*D[i]; FFT(E,-1),FFT(F,-1),FFT(G,-1); for(int i=0;i<nn;i++) ans[i]=((ll)(G[i].real()+0.5)+(((ll)(F[i].real()+0.5)%mod)<<15)+((ll)(E[i].real()+0.5)%mod<<30))%mod; for(ll i=0;i<k;i++) { ans[k*2-i]=(ans[k*2-i]*qpow(k,mod-2)%mod)*qpow(w,c1(i))%mod; printf("%lld\n",ans[k*2-i]); } }