LOJ#3058. 「HNOI2019」白兔之舞

LOJ#3058. 「HNOI2019」白兔之舞

https://loj.ac/problem/3058

分析

• 不妨设转移矩阵为(f)，且用(f^i)来表示走(i)步从(x)到(y)的方案数。

• 跳了(i)步的方案数为(sumlimits_{j=i}^Linom{j-1}{i-1}=inom{L}{i})

• 直接推式子(ans_t=sumlimits_{i=0}^L[imod k=t]f^iinom{L}{i})

• 发现这是个挺套路的单位根反演

• (ans_t imes k=sumlimits_{i=0}^{L}sumlimits_{j=0}^{k-1}w^{j(i-t)}f^iinom{L}{i})

  • ​ $ =sumlimits_{j=0}^{k-1}w^{-jt}sumlimits_{i=0}^Lw^{ij}f^iinom{L}{i}$
  • ​ (=sumlimits_{j=0}^{k-1}w^{-jt}(w^jf+I))

• 把乘积化成和的形式，有(ij=inom{i+j}{2}-inom{i}{2}-inom{j}{2})

• (ans_t imes k=sumlimits_{j=0}^{k-1}w^{-inom{j+t}{2}+inom{j}{2}+inom{t}{2}}g_j)

• (ans_t imes k imes w^{-inom{t}{2}}=sumlimits_{j=0}^{k-1}w^{-inom{j+t}{2}}w^{inom{j}{2}}g_j​)

• 是一个卷积的形式，直接用(mtt)做就行了，用https://cmxrynp.github.io/2019/01/07/fft-optimization/中的(4)次或(3.5)次dft再加上一点点卡常即可通过此题。

代码

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <vector>
#include <cmath>
#include <iostream>
using namespace std;
typedef long long ll;
typedef long double f2;
#define N 262190
#define db(x) cerr<<#x<<" = "<<x<<endl
namespace groot {
	int pri[100],cnt;
	ll qp(ll x,ll y,ll p) {
		ll re=1;for(;y;y>>=1,x=x*x%p)if(y&1)re=re*x%p; return re;
	}
	bool check(int x,int phi) {
		int i;
		for(i=1;i<=cnt;i++) if(qp(x,phi/pri[i],phi+1)==1) return 0;
		return 1;
	}
	int get_g(int p,int phi) {
		int x=phi,i;
		for(i=2;i*i<=x;i++) {
			if(x%i==0) {
				pri[++cnt]=i;
				while(x%i==0) x/=i;
			}
		}
		if(x!=1) pri[++cnt]=x;
		int g=2;while(!check(g,phi)) g++;
		return g;
	}
}
int mod,n,K,L,X,Y,G;
ll g[N],w[N],f[N];
ll qp(ll x,ll y) {
	ll re=1; for(;y;y>>=1,x=x*x%mod) if(y&1) re=re*x%mod; return re;
}
ll C2(ll x) {return x*(x-1)/2%K;}
namespace matrix {
	struct Mat {
		ll a[3][3];
		Mat() {memset(a,0,sizeof(a));}
		Mat operator * (const Mat &u) const {
			int i,j,k; Mat re;
			for(i=0;i<n;i++)for(k=0;k<n;k++)for(j=0;j<n;j++) {
				re.a[i][j]+=a[i][k]*u.a[k][j];
			}for(i=0;i<n;i++)for(j=0;j<n;j++)re.a[i][j]%=mod;
			return re;
		}
		Mat operator * (const ll &x) const {
			Mat re=*this; int i,j;
			for(i=0;i<n;i++)for(j=0;j<n;j++)re.a[i][j]=re.a[i][j]*x%mod;
			return re;
		}
		Mat operator + (const Mat &u) const {
			Mat re=*this; int i,j;
			for(i=0;i<n;i++)for(j=0;j<n;j++)re.a[i][j]=(re.a[i][j]+u.a[i][j])%mod;
			return re;
		}
	}I,O;
	Mat qp(Mat x,int y) {
		Mat re=I;
		for(;y;y>>=1,x=x*x)if(y&1)re=re*x; return re;
	}
	void solve() {
		int i;
		for(i=0;i<n;i++) I.a[i][i]=1;
		for(i=0;i<K;i++) {
			g[i]=qp(O*w[i]+I,L).a[X-1][Y-1]*w[C2(i)]%mod;
		}
	}
}
ll tmp[N];
namespace fft {
	const f2 pi=acos(-1);
	struct cp {
		f2 x,y;
		cp() {}
		cp(f2 x_,f2 y_) {x=x_,y=y_;}
		cp operator + (const cp &u) const {return cp(x+u.x,y+u.y);}
		cp operator - (const cp &u) const {return cp(x-u.x,y-u.y);}
		cp operator * (const cp &u) const {return cp(x*u.x-y*u.y,x*u.y+y*u.x);}
	}A[N],B[N],C[N],D[N];
	void dft(cp *a,int len) {
		int i,j,k,t; cp tmp,w,wn;
		for(i=k=0;i<len;i++) {
			if(i>k) swap(a[i],a[k]);
			for(j=len>>1;(k^=j)<j;j>>=1);
		}
		for(k=2;k<=len;k<<=1) {
			wn=cp(cos(2*pi/k),sin(2*pi/k));
			for(t=k>>1,i=0;i<len;i+=k) {
				for(w=cp(1,0),j=i;j<i+t;j++) {
					tmp=a[j+t]*w; a[j+t]=a[j]-tmp; a[j]=a[j]+tmp; w=w*wn;
				}
			}
		}
	}
	int solve() {
		int i;
		w[K]=1;
		for(i=0;i<K+K;i++) f[i]=w[K-C2(i)];
		reverse(g,g+K+1);
		int n=K+K-1,m=K;
		int l=1;
		while(l<=(n+m))l<<=1;
		//for(i=0;i<l;i++) printf("%lld %lld
",f[i],g[i]);
		for(i=0;i<=n;i++) {
			A[i].x=f[i]>>15;
			A[i].y=f[i]&32767;
		}
		for(i=0;i<=m;i++) {
			B[i].x=g[i]>>15;
			B[i].y=g[i]&32767;
		}
		dft(A,l),dft(B,l);
		for(i=0;i<l;i++) {
			int pl=(l-1)&(l-i);
			const cp ca=cp(A[pl].x,-A[pl].y),cb=cp(B[pl].x,-B[pl].y);
			const cp a=(A[i]+ca)*cp(0.5,0),b=(A[i]-ca)*cp(0,-0.5),c=(B[i]+cb)*cp(0.5,0),d=(B[i]-cb)*cp(0,-0.5);
			C[pl]=a*c+a*d*cp(0,1);
			D[pl]=b*c+b*d*cp(0,1);
		}
		dft(C,l),dft(D,l);
		ll invk=qp(K,mod-2);
		for(i=0;i<l;i++) C[i].x/=l,C[i].y/=l,D[i].x/=l,D[i].y/=l;
		//for(i=0;i<=n;i++) for(int j=0;j<=m;j++) {
			//tmp[i+j]=(tmp[i+j]+f[i]*g[j])%mod;
		//}
		for(i=K;i<K+K;i++) {
			ll a=C[i].x+0.5,b=C[i].y+0.5,c=D[i].x+0.5,d=D[i].y+0.5;
			a=((a%mod)<<30)+(((b+c)%mod)<<15)+d;
			a=(a%mod+mod)%mod;
			//a=tmp[i];
			//db(a);
			ll t=w[C2(i-K)]*invk%mod;
			printf("%lld
",a*t%mod);
		}
		return 0;
	}
}
int main() {
	scanf("%d%d%d%d%d%d",&n,&K,&L,&X,&Y,&mod);
	G=groot::get_g(mod,mod-1);
	w[0]=1;w[1]=qp(G,(mod-1)/K);
	int i,j;
	for(i=2;i<K;i++) w[i]=w[i-1]*w[1]%mod;
	//printf("%d
",G);
	for(i=0;i<n;i++)for(j=0;j<n;j++)scanf("%lld",&matrix::O.a[i][j]);
	matrix::solve();
	return fft::solve();
}