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\) 的方案数(到达结点不同算相同方案)。那么:
考虑到 \(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}\),转移矩阵为:
则有:
本质上,\(g\) 的转移完全没有考虑“层”的跨越,而 \(f\) 则需要考虑。具体地,对于 \(g(i,j)\) 所计数的一种从 \((0,x)\) 出发,\(j\) 步后走向某个 \(\in S_j\) 的点的方案,我们需要为其中每个结点安排一个实际所属的层。即:
记要求的答案为 \(\{q_0,q_1,\cdots,q_{k-1}\}\),那么:
单位根反演 \([i\bmod k]\):
似乎不能化简了,带入 \(f(i,j)=\binom{L}ig(i,j)\),进而带入 \(G(i)\)。令 \(Q_r\) 为 \(1\times n\) 的向量,第 \(t\) 维表示 \(y=t\) 时的 \(q_r\),有:
不出意料地发现 \(\omega_k^{ji}\) 能够被收进前面的幂中,继而可以收拢一个二项式定理形式的和式:
后面的矩阵幂可以暴力求,令 \(h_j\) 为向量 \(G(0)(\omega_k^jS+1)^L\) 的第 \(y\) 维,则答案为:
不能暴力 \(\mathcal O(k^2)\) 求,这里用一种科幻的方法:令 \(\mathcal Q(x)=\sum_{i=0}^{k-1}q_ix^i\),我们希望直接求出 \(\mathcal Q(x)\) 的点值表示,然后给它 IDFT 回去直接求出所有 \(q\)。求点值表示理论上要用卷积,所以得把 \(q_r\) 写称卷积的样子。接下来是 trick 式的变形:
令 \(\mathrm f_i=h_j\omega_k^{\binom{j}2}\),\(\mathrm g_i=\omega_k^{-\binom{r+j}2}\),则:
把 \(\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;
}