loj#3058 「HNOI2019」白兔之舞 MTT+機關根反演+矩陣快速幂

Description

自己看吧

Solution

考場上隻想到了n=1的O(L)做法，死于姿勢不足不會機關根反演.jpg

正解是真的難寫，套一個MTT是最惡心的。。

考慮20分怎麼做，一個比較naive的dp就是f[i,j]表示走了i步到了第j列的方案數，nL2轉移就可以成功簽到了

考慮n=1的情況怎麼做。當w=1時，可以發現走恰好d步的方案就是一個組合數，那麼答案實際上就是 f ( d ) = ∑ i = 0 ⌊ L k ⌋ ( L k i + d ) f\left(d\right)=\sum_{i=0}^{\lfloor\frac{L}{k}\rfloor}{\binom{L}{ki+d}} f(d)=i=0∑⌊kL​⌋​(ki+dL​)

這樣就可以獲得0分的好成績

可以發現上柿等價于 ∑ i = 0 L [ k ∣ ( i − d ) ] ( L i ) w i \sum_{i=0}^{L}{{\left[k|\left(i-d\right)\right]}{\binom{L}{i}}w^{i}} i=0∑L​[k∣(i−d)](iL​)wi

套一個機關根反演就有

1 k ∑ j = 0 k − 1 ω k − d j ∑ i = 0 L W i ω k i j ( L i ) \frac{1}{k}\sum_{j=0}^{k-1}{\omega_k^{-dj}\sum_{i=0}^{L}{W^i{\omega_k}^{ij}\binom{L}{i}}} k1​j=0∑k−1​ωk−dj​i=0∑L​Wiωk​ij(iL​)

可以發現後面那一坨可以二項式定理搞成矩陣快速幂，那麼就隻和j有關了。記 C i = ( ω k j W + I ) L C_i={\left({\omega_k}^{j}W+I\right)}^{L} Ci​=(ωk​jW+I)L

很顯然這裡I是機關矩陣，W是給出的鄰接矩陣。實際操作我們隻需要取出W[x][y]就可以了

現在考慮怎麼算 ∑ i = 0 k − 1 w k − d i C i \sum_{i=0}^{k-1}{{w_k}^{-di}C_i} i=0∑k−1​wk​−diCi​

有個騷操作就是構造 a b = ( a + b 2 ) − ( a 2 ) − ( b 2 ) ab=\binom{a+b}{2}-\binom{a}{2}-\binom{b}{2} ab=(2a+b​)−(2a​)−(2b​)證明的話可以考慮拆開來直接算

那麼上柿子就可以變成

ω k ( d 2 ) ∑ i = 0 k − 1 ω k ( i 2 ) C i ⋅ ω k − ( d + i 2 ) {\omega_k}^{\binom{d}{2}}\sum_{i=0}^{k-1}{\omega_k}^{\binom{i}{2}}C_i\cdot {\omega_k}^{-\binom{d+i}{2}} ωk​(2d​)i=0∑k−1​ωk​(2i​)Ci​⋅ωk​−(2d+i​)

可以發現這樣就能構造卷積了，把其中一個翻轉然後結果的第k+i位對應原本的第i位

Code

#include <stdio.h>
#include <string.h>
#include <algorithm>
#include <math.h>
#define rep(i,st,ed) for (int i=st;i<=ed;++i)
#define fill(x,t) memset(x,t,sizeof(x))

typedef long long LL;
typedef long double ld;
const ld pi=acos(-1);
const int N=524588;
const int M=32767;

struct Mat {
	LL r[4][4];
	Mat() {fill(r,0);}
	LL* operator [](int x) {
		return r[x];
	}
} I,W;

struct com {
	ld r,i;
} ;

int rv[N],v[N],w[N],MOD,n,k;
int A[N],B[N],C[N],D[N];

void upd(LL &x,LL v) {
	x+=v; (x>=MOD)?(x-=MOD):0;
}

com operator +(const com &a,const com &b) {return (com) {a.r+b.r,a.i+b.i};}
com operator -(const com &a,const com &b) {return (com) {a.r-b.r,a.i-b.i};}
com operator *(const com &a,const 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};}

int read() {
	int x=0,v=1; char ch=getchar();
	for (;ch<'0'||ch>'9';v=(ch=='-')?(-1):v,ch=getchar());
	for (;ch<='9'&&ch>='0';x=x*10+ch-'0',ch=getchar());
	return x*v;
}

Mat operator *(Mat A,Mat B) {
	Mat C; rep(i,1,3) rep(k,1,3) {
		rep(j,1,3) upd(C[i][j],A[i][k]*B[k][j]%MOD);
	}
	return C;
}

Mat operator *(Mat A,LL b) {
	Mat C;
	rep(i,1,3) rep(j,1,3) C[i][j]=A[i][j]*b%MOD;
	return C;
}

Mat operator +(Mat A,Mat B) {
	Mat C;
	rep(i,1,3) rep(j,1,3) C[i][j]=(A[i][j]+B[i][j])%MOD;
	return C;
}

Mat ksm(Mat A,int dep) {
	Mat C=I;
	for (;dep;dep>>=1,A=A*A) {
		if (dep&1) C=C*A;
	}
	return C;
}

int ksm(int x,int dep) {
	int res=1;
	for (;dep;dep>>=1,x=(LL)x*x%MOD) {
		(dep&1)?(res=(LL)res*x%MOD):0;
	}
	return res;
}

int P(int n) {
	return ((LL)n*(n-1)/2)%k;
}

void FFT(com *a,int n) {
	for (int i=0;i<n;++i) if (i<rv[i]) std:: swap(a[i],a[rv[i]]);
	for (int i=1;i<n;i<<=1) {
		com wn=(com) {cos(pi/i),sin(pi/i)};
		for (int j=0;j<n;j+=(i<<1)) {
			com w=(com) {1,0};
			for (int k=0;k<i;++k) {
				com u=a[j+k],v=a[j+k+i]*w;
				a[j+k]=u+v,a[j+k+i]=u-v;
				w=w*wn;
			}
		}
	}
}

void mul(int *A,int *B,int *C,int n) {
	static com a[N],b[N],Da[N],Db[N],Dc[N],Dd[N];
	for (int i=0;i<n;++i) A[i]=(A[i]+MOD)%MOD;
	for (int i=0;i<n;++i) B[i]=(B[i]+MOD)%MOD;
	for (int i=0;i<n;++i) a[i]=(com) {(ld)(A[i]&M),(ld)(A[i]>>15)};
	for (int i=0;i<n;++i) b[i]=(com) {(ld)(B[i]&M),(ld)(B[i]>>15)};
	FFT(a,n),FFT(b,n);
	for (int i=0;i<n;++i) {
		int j=(n-1)&(n-i);
		com da=(a[i]+conj(a[j]))*(com) {0.5,0};
		com db=(a[i]-conj(a[j]))*(com) {0,-0.5};
		com dc=(b[i]+conj(b[j]))*(com) {0.5,0};
		com dd=(b[i]-conj(b[j]))*(com) {0,-0.5};
		Da[j]=da*dc,Db[j]=da*dd,Dc[j]=db*dc,Dd[j]=db*dd;
	}
	for (int i=0;i<n;++i) a[i]=Da[i]+Db[i]*(com) {0,1};
	for (int i=0;i<n;++i) b[i]=Dc[i]+Dd[i]*(com) {0,1};
	FFT(a,n),FFT(b,n);
	for (int i=0;i<n;++i) {
		LL da=(LL)(a[i].r/n+0.5)%MOD;
		LL db=(LL)(a[i].i/n+0.5)%MOD;
		LL dc=(LL)(b[i].r/n+0.5)%MOD;
		LL dd=(LL)(b[i].i/n+0.5)%MOD;
		C[i]=((LL)da+((LL)((db+dc)%MOD)<<15)+((LL)(dd%MOD)<<30)%MOD)%MOD;
	}
}

int get_g() {
	int tmp=MOD-1; v[0]=0;
	for (int i=2;i*i<=tmp;++i) {
		if (tmp%i==0) {
			while (tmp%i==0) tmp/=i;
			v[++v[0]]=i;
		}
	}
	if (tmp!=1) v[++v[0]]=tmp;
	rep(g,2,MOD) {
		bool flag=false;
		rep(i,1,v[0]) {
			if (ksm(g,(MOD-1)/v[i])==1) {
				flag=true; break;
			}
		}
		if (!flag) return g;
	}
}

int main(void) {
	freopen("data.in","r",stdin);
	int n=read(); k=read(); int L=read();
	int x=read(),y=read(); MOD=read();
	rep(i,0,n) I[i][i]=1;
	rep(i,1,n) rep(j,1,n) W[i][j]=read();
	int g=get_g();
	w[0]=1,w[1]=ksm(g,(MOD-1)/k);
	rep(i,2,k-1) w[i]=(LL)w[i-1]*w[1]%MOD;
	rep(i,0,k-1) {
		C[i]=ksm(W*w[i]+I,L)[x][y];
		A[i]=(LL)w[P(i)]*C[i]%MOD;
	}
	rep(i,0,k*2-1) B[i]=w[(k-P(i))%k];
	rep(i,0,k/2) std:: swap(A[i],A[k-i]);
	int len=1,lg=0; for (;len<=k*3;) len<<=1,lg++;
	for (int i=0;i<len;++i) rv[i]=(rv[i>>1]>>1)|((i&1)<<(lg-1));
	mul(A,B,D,len);
	rep(i,0,k-1) {
		LL tmp=(LL)ksm(k,MOD-2)*w[P(i)]%MOD;
		LL ans=(LL)D[i+k]*tmp%MOD;
		ans=(ans%MOD+MOD)%MOD;
		printf("%lld\n", ans);
	}
	return 0;
}