【LOJ3058】【HNOI2019】白兔之舞

Description

https://loj.ac/problem/3058

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∑L​Ax,yi​(iL​)k1​j=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 k1​j=0∑k−1​ωk−tj​i=0∑L​(iL​)ωkij​Ax,yi​=k1​j=0∑k−1​ωk−tj​(ωkj​A+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​)​(ωkj​A+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]);
	}
}