(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;
}