#3058. 「HNOI2019」白兔之舞（单位根反演）

题目描述

https://loj.ac/problem/3058

单位根反演

因为ω太难写了所以用w代替

有([n|k]=frac{1}{n}sum_{i=0}^{n-1} w_n^{ik})

证明：

当n|k时显然是1，否则(frac{w_n^{nk}-1}{w^n-1}=0)

题解

一开始想矩乘存多项式然后快速幂循环卷积，然后多乘了一个log直接炸飞

先暴力找p-1的原根，然后即可求出k次单位根w，A是读入的矩阵

(ans_t=sum_{i=0}^L [k|(i-t)]A_{x,y}^i inom{L}{i})

(=frac{1}{k}sum_{i=0}^L sum_{j=0}^{k-1}w^{j(i-t)} A_{x,y}^i inom{L}{i})

(=frac{1}{k}sum_{j=0}^{k-1}w^{-jt} sum_{i=0}^L w^{ij} A_{x,y}^i inom{L}{i})

(=frac{1}{k}sum_{j=0}^{k-1}w^{-jt} (w^jA+1)^L_{x,y})

(=frac{1}{k}sum_{j=0}^{k-1}w^{inom{j}{2}}w^{inom{t}{2}}w^{-inom{j+t}{2}} (w^jA+1)^L_{x,y})

把j变成k-1-j即可卷积，要写mtt

注意mtt求a-bi的点值时的关系是共轭而不是相等

code

#include <bits/stdc++.h>
#define fo(a,b,c) for (a=b; a<=c; a++)
#define fd(a,b,c) for (a=b; a>=c; a--)
#define add(a,b) a=((a)+(b))%mod
#define ll long long
#define D 10000
//#define file
using namespace std;

int N,len,n,K,L,x,y,mod,Mod,i,j,k,l;
ll a[262144],b[262144],g,w,s;

struct type{double x,y;} AA[262144];
type operator + (type a,type b) {return {a.x+b.x,a.y+b.y};}
type operator - (type a,type b) {return {a.x-b.x,a.y-b.y};}
type operator * (type a,type b) {return {a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
struct mat{ll a[3][3]; void clear() {memset(a,0,sizeof(a));}} I,A,B;
mat operator * (mat a,mat b)
{
	static mat c;
	int i,j,k;
	
	fo(i,0,n-1)
	{
		fo(j,0,n-1)
		{
			c.a[i][j]=0;
			fo(k,0,n-1)
			add(c.a[i][j],a.a[i][k]*b.a[k][j]);
		}
	}
	return c;
}
mat Qpower(mat a,ll b) {mat ans=I; while (b) {if (b&1) ans=ans*a;a=a*a;b>>=1;} return ans;}

ll qpower(ll a,ll b) {ll ans=1; while (b) {if (b&1) ans=ans*a%mod;a=a*a%mod;b>>=1;} return ans;}
void dft(type *a,int tp)
{
	int i,j,k,l,S=N,s1=2,s2=1;
	type w,W,u,v;
	
	fo(i,0,N-1)
	{
		j=i,k=0;
		fo(l,1,len)
		k=k*2+(j&1),j>>=1;
		AA[i]=a[k];
	}
	memcpy(a,AA,N*sizeof(type));
	
	fo(i,1,len)
	{
		S>>=1;
		fo(j,0,S-1)
		{
			fo(k,0,s2-1)
			{
				W={cos(2*M_PI*k/s1),sin(2*M_PI*k/s1)*tp};
				u=a[j*s1+k],v=a[j*s1+k+s2]*W;
				a[j*s1+k]=u+v;
				a[j*s1+k+s2]=u-v;
			}
		}
		s1<<=1,s2<<=1;
	}
}
void mul(ll *a,ll *b)
{
	static type A[262144],B[262144],P[262144],Q[262144];
	static ll ans[262144];
	int i,j,k,l;
	
	fo(i,0,N-1) A[i]={a[i]/D,a[i]%D},B[i]={b[i]/D,b[i]%D};
	dft(A,1),dft(B,1);
	fo(i,0,N-1) P[i]=A[i]*B[i],Q[i]=type{A[(N-i)%N].x,-A[(N-i)%N].y}*B[i];
	dft(P,-1),dft(Q,-1);
	fo(i,0,N-1) P[i].x/=N,Q[i].x/=N,P[i].y/=N,Q[i].y/=N;
	
	fo(i,0,N-1)
	ans[i]=(((ll)floor((P[i].x+Q[i].x)/2+0.5)%mod*D%mod*D)+(ll)floor(P[i].y+0.5)*D+(ll)floor((Q[i].x-P[i].x)/2+0.5))%mod;
	memcpy(a,ans,N*8);
}

bool pd(ll g)
{
	int i,j,k,l,s=floor(sqrt(mod-1));
	fo(i,1,s)
	if (!((mod-1)%i))
	{
		if (qpower(g,i)==1) return 0;
		if (i>1 && qpower(g,(mod-1)/i)==1) return 0;
	}
	return 1;
}

int main()
{
	#ifdef file
	freopen("loj3058.in","r",stdin);
	#endif
	
	scanf("%d%d%d%d%d%d",&n,&K,&L,&x,&y,&mod),--x,--y,len=ceil(log2(K*3)),N=pow(2,len),Mod=mod-2;
	I.a[0][0]=I.a[1][1]=I.a[2][2]=1;
	fo(i,0,n-1)
	{
		fo(j,0,n-1)
		scanf("%lld",&A.a[i][j]);
	}
	for (g=1; !pd(g); ++g);
	w=qpower(g,(mod-1)/K);
	
	fo(i,0,K*2-1) a[i]=qpower(qpower(w,1ll*i*(i-1)/2),Mod);
	fo(i,0,K-1)
	{
		B=A,s=qpower(w,i);
		fo(j,0,n-1)
		{
			fo(k,0,n-1)
			B.a[j][k]=B.a[j][k]*s%mod+(j==k);
		}
		B=Qpower(B,L);
		b[(K-1)-i]=qpower(w,1ll*i*(i-1)/2)*B.a[x][y]%mod;
	}
	mul(a,b);
	fo(i,0,K-1)
	printf("%lld
",(a[i+(K-1)]*qpower(K,Mod)%mod*qpower(w,1ll*i*(i-1)/2)%mod+mod)%mod);
	
	fclose(stdin);
	fclose(stdout);
	return 0;
}