Solution -「HNOI2019」「洛谷 P5293」白兔之舞

(mathcal{Description})

  Link.

  不想概括题意.jpg

(mathcal{Solution})

  定义点集 (S_c={(u,v)|v=c})；第 (k) 层点表示所有满足 (u=k) 的结点 ((u,v))。

  尝试朴素 DP，令 (f(i,j)) 表示兔子从 ((0,x)) 出发跳 (i) 步到达某个 ((u,v)in S_j) 的方案数（到达结点不同算不同方案）；(g(i,j)) 表示兔子从 ((0,x)) 出发跳 (i) 步到达任意 ((u,v)in S_j) 的方案数（到达结点不同算相同方案）。那么：

[g(i,j)=sum_{t=0}^mw(t,j)g(i-1,t) ]

  考虑到 (nle3)，令 (G(i)=egin{pmatrix}g(i,0)&g(i,1)&g(i,2)end{pmatrix}) ，初始矩阵 (G(0)=egin{pmatrix}1&0&0end{pmatrix})，转移矩阵为：

[S=egin{pmatrix} w_{11}&cdots&w_{1n}\ vdots&ddots&vdots\ w_{n1}&cdots&w_{nn} end{pmatrix} ]

  则有：

[G(i)=G(0)S^i ]

  本质上，(g) 的转移完全没有考虑“层”的跨越，而 (f) 则需要考虑。具体地，对于 (g(i,j)) 所计数的一种从 ((0,x)) 出发，(j) 步后走向某个 (in S_j) 的点的方案，我们需要为其中每个结点安排一个实际所属的层。即：

[f(i,j)=inom{L}ig(i,j) ]

  记要求的答案为 ({q_0,q_1,cdots,q_{k-1}})，那么：

[q_r=sum_{i=0}^L[imod k=r]f(i,y) ]

  单位根反演 ([imod k])：

[egin{aligned} q_r&=frac{1}ksum_{i=0}^Lf(i,y)sum_{j=0}^{k-1}omega_k^{(i-r)j}\ &=frac{1}ksum_{j=0}^{k-1}omega_k^{-rj}sum_{i=0}^Lf(i,y)omega_k^{ji} end{aligned} ]

  似乎不能化简了，带入 (f(i,j)=inom{L}ig(i,j))，进而带入 (G(i))。令 (Q_r) 为 (1 imes n) 的向量，第 (t) 维表示 (y=t) 时的 (q_r)，有：

[egin{aligned} &q_r=frac{1}ksum_{j=0}^{k-1}omega_k^{-rj}sum_{i=0}^Linom{L}{i}g(i,y)omega^{ji}\ Rightarrow~~~~&Q_k=frac{1}ksum_{j=0}^{k-1}omega_k^{-rj}sum_{i=0}^Linom{L}{i}G(0)S^iomega^{ji} end{aligned} ]

  不出意料地发现 (omega_k^{ji}) 能够被收进前面的幂中，继而可以收拢一个二项式定理形式的和式：

[egin{aligned} Q_k&=frac{1}ksum_{j=0}^{k-1}omega_k^{-rj}G(0)sum_{i=0}^Linom{L}{i}(omega_k^j S)^i\ &=frac{1}ksum_{j=0}^{k-1}omega_k^{-rj}G(0)(omega_k^j S+1)^L end{aligned} ]

  后面的矩阵幂可以暴力求，令 (h_j) 为向量 (G(0)(omega_k^jS+1)^L) 的第 (y) 维，则答案为：

[q_r=frac{1}ksum_{j=0}^{k-1}omega_k^{-rj}h_j ]

  不能暴力 (mathcal O(k^2)) 求，这里用一种科幻的方法：令 (mathcal Q(x)=sum_{i=0}^{k-1}q_ix^i)，我们希望直接求出 (mathcal Q(x)) 的点值表示，然后给它 IDFT 回去直接求出所有 (q)。求点值表示理论上要用卷积，所以得把 (q_r) 写称卷积的样子。接下来是 trick 式的变形：

[egin{aligned} q_r&=frac{1}ksum_{j=0}^{k-1}omega^{-rj}h_j\ &=frac{1}ksum_{j=0}^{k-1}omega_k^{-inom{r+j}2+inom{r}2+inom{j}2}h_j\ &=frac{omega_k^{inom{r}2}}ksum_{j=0}^{k-1}left( h_jomega_k^{inom{j}2} ight)omega_{k}^{-inom{r+j}2} end{aligned} ]

  令 (mathrm f_i=h_jomega_k^{inom{j}2})，(mathrm g_i=omega_k^{-inom{r+j}2})，则：

[q_r=frac{omega_k^{inom{r}2}}ksum_{j=0}^{k-1}mathrm f_jmathrm g_{r+j} ]

  把 (mathrm g) 翻转一下，和式内乘积的下标和就是定值啦！MTT 一下就能求出所有答案。

  综上，瓶颈是求 (k) 个 (3 imes3) 矩阵的 (L) 次幂，复杂度 (mathcal O(3^2klog L))。

(mathcal{Code})

/* Clearink */

#include <cmath>
#include <cstdio>
#include <vector>
#include <cassert>
#include <iostream>
#include <algorithm>

#define rep( i, l, r ) for ( int i = l, rpbound##i = r; i <= rpbound##i; ++i )
#define per( i, r, l ) for ( int i = r, rpbound##i = l; i >= rpbound##i; --i )

typedef long long LL;

inline void wint ( const int x ) {
	if ( 9 < x ) wint ( x / 10 );
	putchar ( x % 10 ^ '0' );
}

const int MAXLEN = 1 << 18, MAXN = 3, MAXK = 65536;
const double PI = acos ( -1 );
int N, K, L, X, Y, M, g, w[MAXK + 5];
int F[MAXLEN + 5], G[MAXLEN + 5], R[MAXLEN + 5];

inline int mul ( const long long a, const int b ) { return a * b % M; }
inline int add ( int a, const int b ) { return ( a += b ) < M ? a : a - M; }
inline int comb2 ( const int a ) { return a < 2 ? 0 : ( a * ( a - 1ll ) >> 1 ) % K; }
inline int mpow ( int a, LL b ) {
	int ret = 1;
	for ( b %= M - 1; b; a = mul ( a, a ), b >>= 1 ) ret = mul ( ret, b & 1 ? a : 1 );
	return ret;
}

namespace MTT {

int rev[MAXLEN + 5];

struct Complex {
	double x, y;
	Complex () {}
	Complex ( const double tx, const double ty ): x ( tx ), y ( ty ) {}
	inline Complex operator + ( const Complex t ) const {
		return Complex ( x + t.x, y + t.y );
	}
	inline Complex operator - ( const Complex t ) const {
		return Complex ( x - t.x, y - t.y );
	}
	inline Complex operator * ( const Complex t ) const {
		return Complex ( x * t.x - y * t.y, x * t.y + y * t.x );
	}
	inline Complex operator / ( const double t ) const {
		return Complex ( x / t, y / t );
	}
} omega[MAXLEN + 5], P[MAXLEN + 5], Q[MAXLEN + 5], C[MAXLEN + 5], D[MAXLEN + 5], E[MAXLEN + 5], F[MAXLEN + 5];

inline void FFT ( const int n, Complex* A, const int tp ) {
	rep ( i, 0, n - 1 ) if ( i < rev[i] ) std::swap ( A[i], A[rev[i]] );
	for ( int i = 2, stp = 1; i <= n; i <<= 1, stp <<= 1 ) {
		for ( int j = 0; j < n; j += i ) {
			rep ( k, 0, stp - 1 ) {
				Complex w ( omega[n / stp * k].x, tp * omega[n / stp * k].y );
				Complex ev ( A[j + k] ), ov ( w * A[j + k + stp] );
				A[j + k] = ev + ov, A[j + k + stp] = ev - ov;
			}
		}
	}
	if ( !~tp ) rep ( i, 0, n - 1 ) A[i] = A[i] / n;
}

inline void initFFT ( const int lg ) {
	int n = 1 << lg;
	rep ( i, 0, n - 1 ) rev[i] = ( rev[i >> 1] >> 1 ) | ( ( i & 1 ) << lg >> 1 );
	for ( int i = 1; i < n; i <<= 1 ) {
		rep ( k, 0, i - 1 ) {
			omega[n / i * k] = Complex ( cos ( PI * k / i ), sin ( PI * k / i ) );
		}
	}
}

inline void polyConv ( const int n, const int m, const int* A, const int* B, int* R ) {
	for ( int i = 0; i < n; ++ i ) P[i] = Complex ( A[i] & 0x7fff, A[i] >> 15 );
	for ( int i = 0; i < m; ++ i ) Q[i] = Complex ( B[i] & 0x7fff, B[i] >> 15 );
	int lg = 0, len = 1;
	for ( ; len < n + m - 1; len <<= 1, ++ lg );
	initFFT ( lg );
	FFT ( len, P, 1 ), FFT ( len, Q, 1 );
	rep ( i, 0, len - 1 ) {
		Complex t ( P[( len - i ) % len].x, -P[( len - i ) % len].y );
		C[i] = ( P[i] + t ) / 2, D[i] = Complex ( 0, 1 ) * ( t - P[i] ) / 2;
	}
	rep ( i, 0, len - 1 ) {
		Complex t ( Q[( len - i ) % len].x, -Q[( len - i ) % len].y );
		E[i] = ( Q[i] + t ) / 2, F[i] = Complex ( 0, 1 ) * ( t - Q[i] ) / 2;
	}
	rep ( i, 0, len - 1 ) {
		Complex c ( C[i] ), d ( D[i] ), e ( E[i] ), f ( F[i] );
		C[i] = c * e, D[i] = c * f + d * e, E[i] = d * f;
		P[i] = C[i] + Complex ( 0, 1 ) * E[i];
	}
	FFT ( len, D, -1 ), FFT ( len, P, -1 );
	rep ( i, 0, n + m - 2 ) {
		int c = ( LL ( P[i].x + 0.5 ) % M + M ) % M,
			d = ( LL ( D[i].x + 0.5 ) % M + M ) % M,
			e = ( LL ( P[i].y + 0.5 ) % M + M ) % M;
		R[i] = ( c + ( 1ll << 15 ) % M * d % M + ( 1ll << 30 ) % M * e % M ) % M;
	}
}

} // namesapce MTT.

struct Matrix {
	int mat[3][3];
	Matrix (): mat {} {}
	inline int* operator [] ( const int k ) { return mat[k]; }
	friend inline Matrix operator * ( Matrix& A, Matrix& B ) {
		Matrix ret;
		rep ( i, 0, 2 ) rep ( k, 0, 2 ) rep ( j, 0, 2 ) {
			ret[i][j] = add ( ret[i][j], mul ( A[i][k], B[k][j] ) );
		}
		return ret;
	}
} S, T;

inline Matrix mpow ( Matrix A, int b ) {
	Matrix ret; ret[0][0] = ret[1][1] = ret[2][2] = 1;
	for ( ; b; A = A * A, b >>= 1 ) if ( b & 1 ) ret = ret * A;
	return ret;
}

inline int getRt () {
	std::vector<LL> fct;
	LL tmp = M - 1;
	for ( int i = 2; 1ll * i * i <= tmp; ++i ) {
		if ( !( tmp % i ) ) {
			fct.push_back ( i );
			for ( ; !( tmp % i ); tmp /= i );
		}
	}
	if ( tmp > 1 ) fct.push_back ( tmp );
	rep ( r, 2, M ) {
		bool flg = true;
		for ( LL f: fct ) {
			if ( !( flg = mpow ( r, ( M - 1 ) / f ) != 1 ) ) {
				break;
			}
		}
		if ( flg ) return r;
	}
	return assert ( false ), -1;
}

inline void make ( Matrix& A, const int w ) {
	rep ( i, 0, N - 1 ) {
		rep ( j, 0, N - 1 ) A[i][j] = mul ( S[i][j], w );
		A[i][i] = add ( A[i][i], 1 );
	}
}

int main () {
	scanf ( "%d %d %d %d %d %d", &N, &K, &L, &X, &Y, &M ), --X, --Y;
	rep ( i, 0, N - 1 ) rep ( j, 0, N - 1 ) scanf ( "%d", &S[i][j] );
	w[0] = 1, w[1] = mpow ( getRt (), ( M - 1 ) / K );
	rep ( i, 2, K - 1 ) w[i] = mul ( w[i - 1], w[1] );
	rep ( i, 0, K - 1 ) {
		make ( T, w[i] ), T = mpow ( T, L );
		F[i] = mul ( T[X][Y], w[comb2 ( i )] );
	}
	rep ( i, 0, K << 1 ) G[i] = w[( K - comb2 ( i ) % K ) % K];
	std::reverse ( G, G + ( K << 1 | 1 ) );
	MTT::polyConv ( K, K << 1 | 1, F, G, R );
	int invk = mpow ( K, M - 2 );
	rep ( i, 0, K - 1 ) {
		int ans = mul ( invk, mul ( R[2 * K - i], w[comb2 ( i )] ) );
		wint ( ans ), putchar ( '
' );
	}
	return 0;
}