Solution -「HNOI 2019」「洛谷 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)$$

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

$$S=\begin{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)=\binom{L}ig(i,j)$$

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

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

  单位根反演 $[i\bmod k]$：

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

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

$$\begin{aligned} &q_r=\frac{1}k\sum_{j=0}^{k-1}\omega_k^{-rj}\sum_{i=0}^L\binom{L}{i}g(i,y)\omega^{ji}\\ \Rightarrow~~~~&Q_k=\frac{1}k\sum_{j=0}^{k-1}\omega_k^{-rj}\sum_{i=0}^L\binom{L}{i}G(0)S^i\omega^{ji} \end{aligned}$$

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

$$\begin{aligned} Q_k&=\frac{1}k\sum_{j=0}^{k-1}\omega_k^{-rj}G(0)\sum_{i=0}^L\binom{L}{i}(\omega_k^j S)^i\\ &=\frac{1}k\sum_{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}k\sum_{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 式的变形：

$$\begin{aligned} q_r&=\frac{1}k\sum_{j=0}^{k-1}\omega^{-rj}h_j\\ &=\frac{1}k\sum_{j=0}^{k-1}\omega_k^{-\binom{r+j}2+\binom{r}2+\binom{j}2}h_j\\ &=\frac{\omega_k^{\binom{r}2}}k\sum_{j=0}^{k-1}\left( h_j\omega_k^{\binom{j}2} \right)\omega_{k}^{-\binom{r+j}2} \end{aligned}$$

  令 $\mathrm f_i=h_j\omega_k^{\binom{j}2}$，$\mathrm g_i=\omega_k^{-\binom{r+j}2}$，则：

$$q_r=\frac{\omega_k^{\binom{r}2}}k\sum_{j=0}^{k-1}\mathrm f_j\mathrm g_{r+j}$$

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

  综上，瓶颈是求 $k$ 个 $3\times3$ 矩阵的 $L$ 次幂，复杂度 $\mathcal O(3^2k\log 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 ( '\n' );
	}
	return 0;
}