题目分析
设 A A A为给定的矩阵。
余数为 t t t时的答案为:
∑ i = 0 L C L i A i [ i   m o d   k = t ] \sum_{i=0}^LC_L^iA^i[i \bmod{k}=t] i=0∑LCLiAi[imodk=t]
已知单位根反演的式子 1 k ∑ i = 0 k − 1 ω k i n = [ k ∣ n ] \frac{1}{k} \sum_{i=0}^{k-1} \omega_k^{in}=[k|n] k1∑i=0k−1ωkin=[k∣n]
操作一下后得到:
1 k ∑ j = 0 k − 1 ( ω k j A + I ) L ω − t j \frac{1}{k} \sum_{j=0}^{k-1} (\omega_k^jA+I)^L \omega^{-tj} k1j=0∑k−1(ωkjA+I)Lω−tj
因为 − t j = C t 2 + C j 2 − C t + j 2 -tj=C_t^2+C_j^2-C_{t+j}^2 −tj=Ct2+Cj2−Ct+j2
所以原式化为
1 k ω k C t 2 ∑ j = 0 k − 1 ( ω k C j 2 ( ω k j A + I ) L ) ( ω k − C t + j 2 ) \frac{1}{k}\omega_k^{C_t^2} \sum_{j=0}^{k-1} (\omega_k^{C_j^2}(\omega_k^jA+I)^L)(\omega_k^{-C_{t+j}^2}) k1ωkCt2j=0∑k−1(ωkCj2(ωkjA+I)L)(ωk−Ct+j2)
设 a k − 1 − i = ω k C i 2 ( ω k i A + I ) L a_{k-1-i}=\omega_k^{C_i^2}(\omega_k^iA+I)^L ak−1−i=ωkCi2(ωkiA+I)L, b i = ω k − C i 2 b_i=\omega_k^{-C_i^2} bi=ωk−Ci2,上面的式子就是个卷积的形式。
用个MTT即可解决。
代码
#include<bits/stdc++.h>
using namespace std;
#define RI register int
int n,K,invK,L,x,y,P,w,g;
struct matrix{int t[4][4];}A,I;
int qm(int x) {return x>=P?x-P:x;}
int ksm(int x,int y) {
int re=1;
for(;y;y>>=1,x=1LL*x*x%P) if(y&1) re=1LL*re*x%P;
return re;
}
inline matrix operator * (matrix A,matrix B) {
matrix C;
for(RI i=0;i<n;++i)
for(RI j=0;j<n;++j) C.t[i][j]=0;
for(RI k=0;k<n;++k)
for(RI i=0;i<n;++i)
for(RI j=0;j<n;++j)
C.t[i][j]=qm(C.t[i][j]+1LL*A.t[i][k]*B.t[k][j]%P);
return C;
}
inline matrix operator + (matrix A,matrix B) {
for(RI i=0;i<n;++i)
for(RI j=0;j<n;++j) A.t[i][j]=qm(A.t[i][j]+B.t[i][j]);
return A;
}
inline matrix operator * (matrix A,int B) {
for(RI i=0;i<n;++i)
for(RI j=0;j<n;++j) A.t[i][j]=1LL*B*A.t[i][j]%P;
return A;
}
matrix ksm(matrix X,int y) {
matrix re=I;
for(;y;y>>=1,X=X*X) if(y&1) re=re*X;
return re;
}
int getg() {
int kP=P-1,js=0,pri[30];
for(RI i=2;i*i<=kP;++i)
if(kP%i==0) {pri[++js]=i;while(kP%i==0) kP/=i;}
if(kP!=1) pri[++js]=kP;
for(RI i=2;;++i) {
int flag=1;
if(ksm(i,P-1)!=1) continue;
for(RI j=1;j<=js;++j)
if(ksm(i,(P-1)/pri[j])==1) {flag=0;break;}
if(flag) return i;
}
}
typedef long double db;
typedef long long LL;
const int N=262150,M=32767;
const db pi=acos(-1);
int ka[N],kb[N],kc[N],rev[N];
struct com{db r,i;}a[N],b[N],Aa[N],Ab[N],Ba[N],Bb[N];
com operator + (com A,com B) {return (com){A.r+B.r,A.i+B.i};}
com operator - (com A,com B) {return (com){A.r-B.r,A.i-B.i};}
com operator * (com A,com B) {return (com){A.r*B.r-A.i*B.i,A.r*B.i+A.i*B.r};}
com conj(com A) {return (com){A.r,-A.i};}
void FFT(com *a,int n) {
for(RI i=0;i<n;++i) if(rev[i]>i) swap(a[i],a[rev[i]]);
for(RI i=1;i<n;i<<=1) {
com wn=(com){cos(pi/i),sin(pi/i)};
for(RI j=0;j<n;j+=(i<<1)) {
com t1,t2,w=(com){1,0};
for(RI k=0;k<i;++k,w=w*wn)
t1=a[j+k],t2=w*a[j+i+k],a[j+k]=t1+t2,a[j+i+k]=t1-t2;
}
}
}
void MTT(int n) {
for(RI i=0;i<n;++i) {
a[i]=(com){(db)(ka[i]&M),(db)(ka[i]>>15)};
b[i]=(com){(db)(kb[i]&M),(db)(kb[i]>>15)};
}
FFT(a,n),FFT(b,n);
for(RI i=0;i<n;++i) {
int j=(n-i)&(n-1);
com kAa=(a[i]+conj(a[j]))*(com){0.5,0};
com kAb=(a[i]-conj(a[j]))*(com){0,-0.5};
com kBa=(b[i]+conj(b[j]))*(com){0.5,0};
com kBb=(b[i]-conj(b[j]))*(com){0,-0.5};
Aa[j]=kAa*kBa,Ab[j]=kAa*kBb,Ba[j]=kAb*kBa,Bb[j]=kAb*kBb;
}
for(RI i=0;i<n;++i)
a[i]=Aa[i]+Ab[i]*(com){0,1},b[i]=Ba[i]+Bb[i]*(com){0,1};
FFT(a,n),FFT(b,n);
for(RI i=0;i<n;++i) {
int kAa=(LL)(a[i].r/n+0.5)%P;
int kAb=(LL)(a[i].i/n+0.5)%P;
int kBa=(LL)(b[i].r/n+0.5)%P;
int kBb=(LL)(b[i].i/n+0.5)%P;
kc[i]=((LL)kAa+((LL)(kAb+kBa)<<15)+((LL)kBb<<30))%P;
}
}
int C2(int x) {return (1LL*x*(x-1)/2)%(P-1);}
int main()
{
scanf("%d%d%d%d%d%d",&n,&K,&L,&x,&y,&P),--x,--y;
for(RI i=0;i<n;++i)
for(RI j=0;j<n;++j) scanf("%d",&A.t[i][j]);
for(RI i=0;i<n;++i) I.t[i][i]=1;
g=getg(),w=ksm(g,(P-1)/K);
for(RI i=0,ww=1;i<K;++i,ww=1LL*ww*w%P) {
matrix tmp=ksm(A*ww+I,L);
ka[K-1-i]=1LL*ksm(w,C2(i))*tmp.t[x][y]%P;
}
for(RI i=0;i<=2*K-2;++i) kb[i]=ksm(ksm(w,C2(i)),P-2);
int kn=1,len=0;
while(kn<=3*K-3) kn<<=1,++len;
for(RI i=0;i<kn;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
MTT(kn),invK=ksm(K,P-2);
for(RI i=0;i<K;++i) printf("%lld\n",1LL*invK*ksm(w,C2(i))%P*kc[K+i-1]%P);
return 0;
}