Description
Solution
首先答案长这样子:
a
n
s
t
=
∑
i
=
0
L
[
k
∣
(
i
−
t
)
]
A
i
(
L
i
)
ans_t=\sum\limits_{i=0}^L[k|(i-t)]A^i\binom{L}{i}
anst=i=0∑L[k∣(i−t)]Ai(iL),
A
A
A是读入的矩阵,最后取
(
x
,
y
)
(x,y)
(x,y)的值
套上单位根反演就是:
∑
i
=
0
L
A
x
,
y
i
(
L
i
)
1
k
∑
j
=
0
k
−
1
ω
k
(
i
−
t
)
j
\sum\limits_{i=0}^LA^i_{x,y}\binom{L}{i}\dfrac{1}{k}\sum\limits_{j=0}^{k-1}\omega_k^{(i-t)j}
i=0∑LAx,yi(iL)k1j=0∑k−1ωk(i−t)j,其中
ω
k
=
g
p
−
1
k
\omega_k=g^{\frac{p-1}{k}}
ωk=gkp−1,
g
g
g是原根。
交换主体再化简:
1
k
∑
j
=
0
k
−
1
ω
k
−
t
j
∑
i
=
0
L
(
L
i
)
ω
k
i
j
A
x
,
y
i
=
1
k
∑
j
=
0
k
−
1
ω
k
−
t
j
(
ω
k
j
A
+
1
)
L
\dfrac{1}{k}\sum\limits_{j=0}^{k-1}\omega_k^{-tj}\sum\limits_{i=0}^L\binom{L}{i}\omega_k^{ij}A_{x,y}^i=\dfrac{1}{k}\sum\limits_{j=0}^{k-1}\omega_k^{-tj}(\omega_{k}^jA+1)^L
k1j=0∑k−1ωk−tji=0∑L(iL)ωkijAx,yi=k1j=0∑k−1ωk−tj(ωkjA+1)L
然后把
−
t
j
-tj
−tj拆开:
−
t
j
=
(
t
2
)
+
(
j
2
)
−
(
t
+
j
2
)
-tj=\binom{t}{2}+\binom{j}{2}-\binom{t+j}{2}
−tj=(2t)+(2j)−(2t+j)。
最后式子就是:
a
n
s
t
=
1
k
ω
k
(
t
2
)
∑
j
=
0
k
−
1
ω
k
(
j
2
)
(
ω
k
j
A
+
1
)
L
ω
k
(
t
+
j
2
)
ans_t=\dfrac{1}{k}\omega_k^{\binom{t}{2}}\sum\limits_{j=0}^{k-1}\omega_k^{\binom{j}{2}}(\omega_k^jA+1)^L\omega_k^{\binom{t+j}{2}}
anst=k1ωk(2t)j=0∑k−1ωk(2j)(ωkjA+1)Lωk(2t+j)
把
ω
k
(
t
+
j
2
)
\omega_k^{\binom{t+j}{2}}
ωk(2t+j)和
a
n
s
t
ans_t
anst翻转一下就是一个卷积形式了。
Code
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#include<cmath>
#define fo(i,j,k) for(int i=j;i<=k;++i)
#define fd(i,j,k) for(int i=j;i>=k;--i)
using namespace std;
typedef long long ll;
const int N=1e5+10,M=3e5+10,qmo=31622;
int n,P;
int w[N];
int qpow(int x,int y){
int s=1;
for(;y;y>>=1,x=(ll)x*x%P) if(y&1) s=(ll)s*x%P;
return s;
}
void inc(int &x,int y){
x=x+y>=P?x+y-P:x+y;
}
struct mat{
int a[4][4];
void clear(){
memset(a,0,sizeof(a));
fo(i,1,n) a[i][i]=1;
}
friend mat operator *(mat x,mat y){
mat z;
memset(z.a,0,sizeof(z.a));//z.a[1] z.a[2] z.a[3]
fo(i,1,n)
fo(j,1,n){
ll t=0;
fo(k,1,n) t+=(ll)x.a[i][k]*y.a[k][j];
z.a[i][j]=t%P;
}
return z;
}
void mul(int x){
fo(i,1,n)
fo(j,1,n) a[i][j]=(ll)a[i][j]*x%P;
}
}A,S,z;
void matpow(int y){
S.clear();
for(;y;y>>=1,A=A*A)
if(y&1) S=S*A;
}
int pr[50];
int getrt(){
int x=P-1;
for(int i=2;i*i<=x;++i) if(x%i==0){
pr[++pr[0]]=i;
for(;x%i==0;x/=i);
}
fo(i,2,P-1){
bool flag=1;
fo(j,1,pr[0]) if(qpow(i,(P-1)/pr[j])==1){
flag=0;
break;
}
if(flag) return i;
}
}
struct node{
double x,y;
node(double _x=0,double _y=0) {x=_x,y=_y;}
}a[M],b[M],c1[M],c2[M],c3[M],W[M];
const double pi=acos(-1);
node operator+(node x,node y) {return node(x.x+y.x,x.y+y.y);}
node operator-(node x,node y) {return node(x.x-y.x,x.y-y.y);}
node operator*(node x,node y) {return node(x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x);}
node conj(node a) {return node(a.x,-a.y);}
int rev[M],g[M];
void DFT(node *a,int fn,int sig) {
fo(i,0,fn-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int m=2;m<=fn;m<<=1){
int half=m/2,t=fn/m;
fo(i,0,half-1){
node w=sig>0?W[t*i]:W[fn-t*i];
for(int j=i;j<fn;j+=m){
node u=a[j],v=w*a[j+half];
a[j]=u+v,a[j+half]=u-v;
}
}
}
if(sig==-1) fo(i,0,fn-1) a[i].x/=fn;
}
void mul(int *A,int *B,int nl,int nr){
int fn,cnt=0;
for(fn=1;fn<=nl+nr;fn<<=1) ++cnt;
fo(i,0,fn) W[i]=node(cos(2*i*pi/fn),sin(2*i*pi/fn));
fo(i,1,fn-1) rev[i]=rev[i>>1]>>1|(i&1)<<(cnt-1);
fo(i,0,nl) a[i]=node(A[i]/qmo,A[i]%qmo);
fo(i,nl+1,fn-1) a[i]=node(0,0);
fo(i,0,nr) b[i]=node(B[i]/qmo,B[i]%qmo);
fo(i,nr+1,fn-1) b[i]=node(0,0);
DFT(a,fn,1),DFT(b,fn,1);
fo(i,0,fn-1){
int j=(fn-i)&(fn-1);
node p1,p2,q1,q2;
p1=(a[i]+conj(a[j]))*node(0.5,0);
p2=(a[i]-conj(a[j]))*node(0,-0.5);
q1=(b[i]+conj(b[j]))*node(0.5,0);
q2=(b[i]-conj(b[j]))*node(0,-0.5);
c1[i]=p1*q1,c2[i]=p1*q2+p2*q1,c3[i]=p2*q2;
}
DFT(c1,fn,-1),DFT(c2,fn,-1),DFT(c3,fn,-1);
fo(i,0,fn-1){
g[i]=((ll)(c1[i].x+0.5)*qmo%P*qmo%P+
(ll)(c2[i].x+0.5)*qmo%P+
(ll)(c3[i].x+0.5))%P;
}
fo(i,0,nl+nr) A[i]=g[i];
fo(i,nl+nr+1,fn-1) A[i]=0;
}
int k;
int C2(int x){
return (ll)x*(x-1)/2%k;
}
int c[M],d[M];
int main()
{
freopen("dance.in","r",stdin);
freopen("dance.out","w",stdout);
int L,x,y;
scanf("%d %d %d %d %d %d",&n,&k,&L,&x,&y,&P);
int nk=qpow(k,P-2);
w[0]=1,w[1]=qpow(getrt(),(P-1)/k);
fo(i,2,k-1) w[i]=(ll)w[i-1]*w[1]%P;
fo(i,1,n)
fo(j,1,n) scanf("%d",&z.a[i][j]);
fo(i,0,k-1){
A=z;
A.mul(w[i]);
fo(j,1,n) inc(A.a[j][j],1);
matpow(L);
c[i]=(ll)S.a[x][y]*w[C2(i)]%P;
}
fo(i,0,2*k-2) d[i]=qpow(w[C2(i)],P-2);
/*fo(i,0,k-1){
int tmp=0;
fo(j,0,k-1) inc(tmp,(ll)c[j]*d[i+j]%P);
tmp=(ll)tmp*nk%P*w[C2(i)]%P;
printf("%d\n",tmp);
}*/
reverse(d,d+2*k-1);
mul(c,d,k-1,2*k-2);
reverse(c,c+2*k-1);
fo(i,0,k-1){
c[i]=(ll)c[i]*nk%P*w[C2(i)]%P;
printf("%d\n",c[i]);
}
}