LGP5293题解

感觉比想象中的简单多了。。。除了比较缝合之外没什么难点。。。

显然有 DP：

设 \(f[t][i][k]\) 表示目前在第 \(t\) 行 \(i\) 列，走的步数对 \(k\) 取模的方案数，\(g[t][i][k]=\sum_{j=0}^{i}f[t][j][k]\)。

显然有 \(f[t][i][k]=\sum_{x=1}^{n}g[x][i-1][k-1\bmod K]\times w[p][q]\)

直接把 \(k\) 的维度看成多项式，显然是一个循环卷积，有：

\[f[q][i](x)=\sum_{p=1}^{n}g[p][i-1](x)\times w[p][q](x) \]

\[g[q][i](x)=g[q][i-1](x)+\sum_{p=1}^{n}g[p][i-1](x)\times w[p][q](x) \]

对于循环卷积最好的办法是 FFT，注意到 \(k\mid(p-1)\) 所以可以使用单位根：

\[g[q][i](\omega_k^t)=g[q][i-1](\omega_k^t)+\sum_{p=1}^{n}g[p][i-1](\omega_k^t)\times w[x][t](\omega_k^t) \]

只需要在最开始的时候对所有相关的东西做一遍长度为 \(k\) 的 FFT 即可，这里使用 Z 变换即可。

继续，每次对 \(\omega_k^t\) 相同的部分做：

\[g_q[i]=g_q[i-1]+\sum_{p=1}^{n}g_p[i-1]\times w[p][q] \]

这里就有常系数齐次线性递推的影子了。

这里稍微推一下可以用常系数齐次线性递推，或者偷懒用矩快也行。

复杂度 \(O(k\log k+kn^3\log L)\)，稳的。

以及这题任意模数需要 MTT，不懂为什么需要单位根反演。

#include<cstdio>
#include<cctype>
#include<cmath>
#define IMP(lim,act) for(int qwq=(lim),i=0;i^qwq;++i)act
typedef double db;
const db Pi=acos(-1);
const int M=1<<16|5;
inline int read(){
	int n(0);char s;while(!isdigit(s=getchar()));while(n=n*10+(s&15),isdigit(s=getchar()));return n;
}
inline void write(int n){
	static char s[15];int top(0);while(s[++top]=n%10^48,n/=10);while(putchar(s[top]),--top);putchar(10);
}
struct Barrett{
	typedef unsigned long long ull;
	typedef __uint128_t LL;
	ull B,m;
	Barrett(const ull&m=2):m(m),B((LL(1)<<64)/m){}
	friend inline ull operator%(const ull&a,const Barrett&mod){
		ull r=a-mod.m*(LL(mod.B)*a>>64);return r>=mod.m?r-mod.m:r;
	}
}mod;
int n,k,L,x,y,g,P,W[3][3],f[M];
int w1,c[M<<1],ic[M<<1],F[M<<1],G[M<<1],H[M<<1];
inline void swap(int&a,int&b){
	int c=a;a=b;b=c;
}
inline int qpow(int a,int b){
	int ans(1);for(;b;b>>=1,a=1ll*a*a%mod)if(b&1)ans=1ll*ans*a%mod;return ans;
}
inline int findg(int P){
	static int p[10];int len(0),phi=P-1;
	for(int i=2;i*i<=phi;++i)if(!(phi%i)){
		p[++len]=i;while(!(phi%i))phi/=i;
	}
	if(phi^1)p[++len]=phi;
	int g=2;
	while(true){
		bool typ(true);
		for(int i=1;i<=len;++i)if(qpow(g,(P-1)/p[i])==1){
			typ=false;break;
		}
		if(typ)return g;++g;
	}
	return-1;
}
struct complex{
	db x,y;
	complex(const db&x=0,const db&y=0):x(x),y(y){}
	inline complex operator+(const complex&it)const{
		return complex(x+it.x,y+it.y);
	}
	inline complex operator-(const complex&it)const{
		return complex(x-it.x,y-it.y);
	}
	inline complex operator*(const complex&it)const{
		return complex(x*it.x-y*it.y,x*it.y+y*it.x);
	}
}buf[M<<2],*w[20];
inline int Getlen(const int&n){
	int len(0);while((1<<len)<n)++len;return len;
}
inline void swap(complex&a,complex&b){
	complex c=a;a=b;b=c;
}
inline int Get(const db&x){
	return((long long)floor(x+.5))%mod;
}
inline void init(const int&n){
	const int&m=Getlen(n);complex*now=buf;w[m]=now;now+=1<<m;
	IMP(1<<m,w[m][i]=complex(std::cos(i*Pi/(1<<m)),std::sin(i*Pi/(1<<m))));
	for(int k=m-1;k>=0&&(w[k]=now,now+=1<<k);--k)IMP(1<<k,w[k][i]=w[k+1][i<<1]);
}
inline void DFT(complex*f,const int&M){
	const int&n=1<<M;
	for(int len=n>>1,d=M-1;d>=0;--d,len>>=1)for(int k=0;k^n;k+=len<<1){
		complex*W=w[d],*L=f+(k),*R=f+(k|len),x,y;IMP(len,(x=*L,y=*R)),*L++=(x+y),*R++=*W++*(x-y);
	}
}
inline void IDFT(complex*f,const int&M){
	const int&n=1<<M;
	for(int len=1,d=0;d^M;++d,len<<=1)for(int k=0;k^n;k+=len<<1){
		complex*W=w[d],*L=f+(k),*R=f+(k|len),x,y;IMP(len,(x=*L,y=*W++**R)),*L++=(x+y),*R++=(x-y);
	}
	IMP(n,(f[i].x/=n,f[i].y/=n));for(int i=1;(i<<1)<n;++i)swap(f[i],f[n-i]);
}
inline void MTT(int*f,int*g,int*h,const int&n,const int&m,const int&LEN){
	static complex Q[M<<1],P[M<<1],T[M<<1];const int&len=Getlen(LEN);
	IMP(n,(Q[i].x=f[i]&32767,P[i].x=f[i]>>15));IMP(m,T[i]=complex(g[i]&32767,g[i]>>15));
	DFT(Q,len);DFT(P,len);DFT(T,len);IMP(1<<len,(Q[i]=Q[i]*T[i],P[i]=P[i]*T[i]));IDFT(Q,len);IDFT(P,len);
	IMP(LEN,h[i]=(Get(Q[i].x)+(1ll*Get(Q[i].y+P[i].x)<<15)+(1ll*Get(P[i].y)<<30))%mod);
	IMP(1<<len,Q[i]=P[i]=T[i]=0);
}
inline void czt_init(const int&n){
	c[0]=ic[0]=1;c[1]=w1=qpow(g=findg(P),(P-1)/n);ic[1]=qpow(c[1],P-2);
	for(int i=2;i^n+n-1;++i)c[i]=1ll*c[i-1]*c[1]%mod;for(int i=2;i^n;++i)ic[i]=1ll*ic[i-1]*ic[1]%mod;
	for(int i=1;i^n+n-1;++i)c[i]=1ll*c[i-1]*c[i]%mod;for(int i=1;i^n;++i)ic[i]=1ll*ic[i-1]*ic[i]%mod;
	IMP(n+n-1,G[i]=c[i]);
}
inline void CZT(int*f,const int&n,const bool&typ){
	IMP(n,F[n-i]=1ll*f[i]*ic[i])%mod;
	MTT(F,G,H,n+1,n+n-1,n+n);IMP(n,f[i]=1ll*H[n+i]*ic[i]%mod);IMP(n+n,F[i]=H[i]=0);
	if(typ)for(int i=1;(i<<1)<n;++i)swap(f[i],f[n-i]);const int&k=qpow(n,P-2);IMP(n,f[i]=1ll*f[i]*k%mod);
}
namespace MATRIX{
	const int S=3;
	struct Matrix{
		int a[S][S];
		Matrix(){
			for(int i=0;i^n;++i)for(int j=0;j^n;++j)a[i][j]=0;
		}
		inline Matrix&I(){
			for(int i=0;i^n;++i)for(int j=0;j^n;++j)a[i][j]=i==j;return*this;
		}
		inline int*operator[](const int&x){
			return a[x];
		}
		inline Matrix operator*(Matrix b){
			static long long c[S][S];Matrix ans;
			for(int i=0;i^n;++i)for(int j=0;j^n;++j)for(int k=0;k^n;++k)c[i][j]+=1ll*a[i][k]*b[k][j];
			for(int i=0;i^n;++i)for(int j=0;j^n;++j)ans[i][j]=c[i][j]%mod,c[i][j]=0;return ans;
		}
		inline Matrix pow(int b){
			Matrix ans,base=*this;ans.I();for(;b;b>>=1,base=base*base)if(b&1)ans=ans*base;return ans;
		}
	};
}
signed main(){
	MATRIX::Matrix m;
	n=read();k=read();L=read();x=read()-1;y=read()-1;mod=Barrett(P=read());init(k<<1);czt_init(k);
	for(int i=0;i^n;++i)for(int j=0;j^n;++j)W[i][j]=read();for(int x=0;x^n;++x)for(int y=0;y^n;++y)m[x][y]=W[x][y];
	for(int i=0;i^k;++i){
		for(int x=0;x^n;++x)++m[x][x];f[i]=m.pow(L)[x][y];for(int x=0;x^n;++x)--m[x][x];
		for(int x=0;x^n;++x)for(int y=0;y^n;++y)m[x][y]=1ll*m[x][y]*w1%mod;
	}
	CZT(f,k,true);IMP(k,write(f[i]));
}