Loj 3058. 「HNOI2019」白兔之舞

Loj 3058. 「HNOI2019」白兔之舞

题目描述

有一张顶点数为 ((L+1) imes n) 的有向图。这张图的每个顶点由一个二元组 ((u,v)) 表示 ((0le ule L,1le vle n))。这张图不是简单图，对于任意两个顶点 ((u_1,v_1),(u_2,v_2))，如果 (u_1<u_2)，则从 ((u_1,v_1)) 到 ((u_2,v_2)) 一共有 (w(v_1,v_2)) 条不同的边，如果 (u_1ge u_2) 则没有边。

白兔将在这张图上上演一支舞曲。白兔初始时位于该有向图的顶点 ((0,x))。

白兔将会跳若干步。每一步，白兔会从当前顶点沿任意一条出边跳到下一个顶点。白兔可以在任意时候停止跳舞（也可以没有跳就直接结束）。当到达第一维为 (L) 的顶点就不得不停止，因为该顶点没有出边。

假设白兔停止时，跳了 (m) 步，白兔会把这只舞曲给记录下来成为一个序列。序列的第 (i) 个元素为它第 (i) 步经过的边。

问题来了：给定正整数 (k) 和 (y (1le yle n))，对于每个 (t (0le t<k))，求有多少种舞曲（假设其长度为 (m)）满足 (m mod k=t)，且白兔最后停在了坐标第二维为 (y) 的顶点？

两支舞曲不同定义为它们的长度（(m)）不同或者存在某一步它们所走的边不同。

输出的结果对 (p) 取模。

输入格式

第一行六个用空格隔开的整数 (n,k,L,x,y,p)。

接下来 (n) 行，每行有 (n) 个用空格隔开的整数，第 (i) 行的第 (j) 个数表示 (w(i,j))。

输出格式

依次输出 (k) 行，每行一个数表示答案对 (p) 取模的结果。

样例

样例输入 1

2 2 3 1 1 998244353
2 1
1 0

样例输出 1

16
18

样例输入 2

3 4 100 1 3 998244353
1 1 1
1 1 1
1 1 1

样例输出 2

506551216
528858062
469849094
818871911

数据范围与提示

对于全部数据，(p) 为一个质数，(10^8<p<2^{30})，(1le nle 3)，(1le xle n)，(1le yle n)，(0le w(i,j)<p)，(1le kle 65536)，(k) 为 (p-1) 的约数，(1le Lle 10^8)。

对于每组测试点，特殊限制如下：

• 测试点 (1,2)：(Lle 10^5)；
• 测试点 (3)：(n=1,w(1,1)=1)，(k) 的最大质因子为 (2)；
• 测试点 (4)：(n=1)，(k) 的最大质因子为 (2)；
• 测试点 (5)：(n=1,w(1,1)=1)；
• 测试点 (6)：(n=1)；
• 测试点 (7,8)：(k) 的最大质因子为 (2)。

(\)

湖南的题好神啊，只能(Orz)了。

先考虑从(x)走(L)步到(y)，假设每次向右只走一格的方案数。很容易写出一个(DP)方程，然后用矩阵快速幂优化。我们读进来的矩阵(A)就是转移矩阵。所以方案数就是([A^i]_{x,y})。

计算答案的时候就枚举走了(i)步：

[ans_t=sum_{i=0}^L[i\%k==t]inom{L}{i}[A^i]_{x,y}\ =sum_{i=0}^Lfrac{1}{k}sum_{j=0}^{k-1}omega_k^{(i-t)*j}inom{L}{i}[A^i]_{x,y}\ =frac{1}{k}sum_{j=0}^{k-1}omega_k^{-jt} sum_{i=0}^Linom{L}{i}(omega_k^{ij} [A^i]_{x,y})\ =frac{1}{k}sum_{j=0}^{k-1}omega_k^{-jt}[(omega_k^jA+I)^L]_{x,y}\ ]

其中：

[I=egin{bmatrix} 1&0&0\ 0&1&0\ 0&0&1\ end{bmatrix} ]

上面最后一步变换实质上就是二项式定理(sum_{i=0}^ninom{n}{i}x_i=(1+x)^n)的运用。

设：

[f_j=[(omega_k^jA+I)^L]_{x,y} ]

则：

[ans_t=frac{1}{k}sum_{j=0}^{k-1}omega_{k}^{-jt}f_j\ ]

然后就是一个很牛逼的套路：

[-tj=inom{t}{2}+inom{j}{2}-inom{t+j}{2} ]

所以：

[ans_t=frac{1}{k}sum_{j=0}^{k-1}omega_k^{inom{t}{2}+inom{j}{2}-inom{t+j}{2}}f_j\ =frac{1}{k}omega_k^{inom{t}{2}}sum_{j=0}^{k-1}omega_k^{inom{j}{2}}f_jcdotomega_k^{-inom{t+j}{2}} ]

(ans_i)的生成函数可以写成下面两个函数卷积。

[F(n)=omega_k^{inom{i}{2}}cdot f_n\ G(n)=omega_k^{-inom{2*k-1-n}{2}}\ ]

代码：

#include<bits/stdc++.h>
#define ll long long
#define K 200005

using namespace std;
inline ll Get() {ll x=0,f=1;char ch=getchar();while(ch<'0'||ch>'9') {if(ch=='-') f=-1;ch=getchar();}while('0'<=ch&&ch<='9') {x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}return x*f;}

ll n,k,L,x,y,mod;
ll g,W,w[K];
ll f[K];
ll ksm(ll t,ll x) {
	ll ans=1;
	for(;x;x>>=1,t=t*t%mod)
		if(x&1) ans=ans*t%mod;
	return ans;
}

vector<int>fac;
void work(int t) {
	for(int i=2,maxx=sqrt(t);i<=maxx;i++) {
		if(t%i==0) fac.push_back(i);
		while(t%i==0) t/=i;
	}
	if(t>1) fac.push_back(t);
}

ll Find_root() {
	for(int i=2;i<mod;i++) {
		if(ksm(i,mod-1)!=1) continue ;
		int flag=1;
		for(int j=0;j<fac.size();j++) {
			if(ksm(i,(mod-1)/fac[j])==1) {
				flag=0;
				break;
			}
		}
		if(flag) return i;
	}
}

struct matrix {
	ll a[4][4];
	void Init() {memset(a,0,sizeof(a));}
}A,B;

matrix operator *(const matrix &a,const matrix &b) {
	static matrix tem;
	tem.Init();
	for(int i=1;i<=n;i++) {
		for(int j=1;j<=n;j++) {
			for(int k=1;k<=n;k++)
				tem.a[i][j]+=a.a[i][k]*b.a[k][j]%mod;
				tem.a[i][j]%=mod;
		}
	}
	return tem;
}

matrix ksm(matrix g,ll x) {
	static matrix ans;
	ans.Init();
	for(int i=1;i<=n;i++) ans.a[i][i]=1;
	for(;x;x>>=1,g=g*g)
		if(x&1) ans=ans*g;
	return ans;
}

ll C_2(ll n) {return n*(n-1)/2%k;}
struct Com {
	long double r,v;
	Com() {}
	Com(long double _r,long double _v) {r=_r,v=_v;}
};

Com operator *(const Com &a,const Com &b) {return Com(a.r*b.r-a.v*b.v,a.r*b.v+a.v*b.r);}
Com operator +(const Com &a,const Com &b) {return Com(a.r+b.r,a.v+b.v);}
Com operator -(const Com &a,const Com &b) {return Com(a.r-b.r,a.v-b.v);}
Com operator /(const Com &a,const long double &b) {return Com(a.r/b,a.v/b);}
const long double pi=acos(-1);
void FFT(Com *a,int d,int flag) {
	int n=1<<d;
	static int rev[K<<2];
	for(int i=0;i<n;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<d-1);
	for(int i=0;i<n;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
	for(int s=1;s<=d;s++) {
		int len=1<<s,mid=len>>1;
		Com w(cos(2*flag*pi/len),sin(2*flag*pi/len));
		for(int i=0;i<n;i+=len) {
			Com t(1,0);
			for(int j=0;j<mid;j++,t=t*w) {
				Com u=a[i+j],v=a[i+j+mid]*t;
				a[i+j]=u+v;
				a[i+j+mid]=u-v;
			}
		}
	}
	if(flag==-1) {for(int i=0;i<n;i++) a[i]=a[i]/n;}
}

ll F[K<<2],G[K<<2],ans[K<<2];
void mul(ll *a,ll *b,ll *c,int d) {
	Com A[K<<2],B[K<<2],C[K<<2],D[K<<2];
	Com f1[K<<2],f2[K<<2],f3[K<<2],f4[K<<2];
	ll maxx=sqrt(mod);
	int n=1<<d;
	for(int i=0;i<n;i++) {
		A[i]=Com(a[i]/maxx,0);
		B[i]=Com(a[i]%maxx,0);
		C[i]=Com(b[i]/maxx,0);
		D[i]=Com(b[i]%maxx,0);
	}
	FFT(A,d,1),FFT(B,d,1);
	FFT(C,d,1),FFT(D,d,1);
	for(int i=0;i<n;i++) {
		f1[i]=A[i]*C[i];
		f2[i]=A[i]*D[i];
		f3[i]=B[i]*C[i];
		f4[i]=B[i]*D[i];
	}
	FFT(f1,d,-1),FFT(f2,d,-1);
	FFT(f3,d,-1),FFT(f4,d,-1);
	for(int i=0;i<n;i++) {
		c[i]=((ll)(f1[i].r+0.5)*maxx%mod*maxx%mod+(ll)(f2[i].r+0.5)*maxx%mod+(ll)(f3[i].r+0.5)*maxx+(ll)(f4[i].r+0.5))%mod;
	}
}

int main() {
	n=Get(),k=Get(),L=Get(),x=Get(),y=Get(),mod=Get();
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++)
			A.a[i][j]=Get();
	
	work(mod-1);
	g=Find_root();
	
	w[0]=1;
	w[1]=ksm(g,(mod-1)/k);
	for(int i=2;i<k;i++) w[i]=w[i-1]*w[1]%mod;
	
	for(int j=0;j<k;j++) {
		B=A;
		for(int i=1;i<=n;i++) {
			for(int k=1;k<=n;k++) {
				B.a[i][k]=B.a[i][k]*w[j];
				if(i==k) B.a[i][k]++;
				B.a[i][k]%=mod;
			}
		}
		B=ksm(B,L);
		f[j]=B.a[x][y];
	}
	
	for(int i=0;i<k;i++) F[i]=w[C_2(i)]*f[i]%mod;
	for(int i=0;i<2*k;i++) G[i]=w[(k-C_2(i))%k]%mod;
	
	reverse(G,G+2*k);
	
	int d=ceil(log2(3*k));
	mul(F,G,ans,d);
	
	ll invk=ksm(k,mod-2);
	for(int i=0;i<k;i++) {
		cout<<ans[2*k-1-i]*invk%mod*w[C_2(i)]%mod<<"
";
	}
	return 0;
}