[HNOI2019] 白兔之舞 题解

白兔虐我千百遍，我待白兔如初恋。

随着单位根反演的普及，这道题或许会逐渐成为模板。

设 \(F_i\) 表示走 \(i\) 步到 \(y\) 得到方案，那么显然有： \[ F_i={L \choose i}A^i[x][y] \] 其中 \(A\) 矩阵是给出的邻接矩阵 \(w\)。（如果这个都不会还学什么单位根反演？）

设 \(G_x\) 表示表示最终的 \(x\) 的答案，那么显然 \(G_x=\sum_iF_i[i\text{ mod }k=x]\)。括号内的内容等价于 \(k|i-x\)。

直接套单位根反演，就有 \(G_x=\frac{1}{k}\sum_{j=0}^{k-1}\omega_k^{-xj}\sum_{i=0}^L{L \choose i}A^i[x][y]\omega_k^{ij}\)。

不难发现后面一堆就是 \((A\omega_k^j+I)^L[x][y]\)，直接矩阵快速幂算即可，设算出来的结果为 \(H_j\)，那么原式就等于 \(\frac{1}{k}\sum_{j=0}^{k-1}\omega_k^{-xj}H_j\)。

按照 Bluestein 的套路，把 \(-xj\) 拆成 \({j \choose 2}+{x \choose 2}-{j+x \choose 2}\)，那么原式就为 \(\frac{1}{k}\omega_k^{x\choose 2}\sum_{j=0}^{k-1}\omega_k^{j \choose 2}H_j\cdot\omega_k^{j+x\choose 2}\)，不难发现右边是一个减法卷积，直接任意模数 FFT 即可。

#include <bits/stdc++.h>
#define RI register int
typedef long long LL;

#define FILEIO(name) freopen(name".in", "r", stdin), freopen(name".out", "w", stdout);

using namespace std;

int const MAXN = 3e5 + 5;
LL n, k, L, x, y, mod;

LL qpow(LL a, LL k) {
  LL re = 1;
  for (; k; k >>= 1, a = a * a % mod)
    if (k & 1) re = re * a % mod;
  return re;
}

int Prime[105];
LL Getwn() {
  int cnt = 0, tmp = mod - 1, sq = sqrt(mod);
  for (RI i = 2; i <= sq; ++i)
    if (tmp % i == 0) {
      Prime[++cnt] = i;
      while (tmp % i == 0)
        tmp /= i;
    }
  if (tmp != 1) Prime[++cnt] = tmp;
  for (RI i = 2; 666; ++i) {
    int flag = 1;
    for (RI j = 1; flag && j <= cnt; ++j)
      if (qpow(i, (mod - 1) / Prime[j]) == 1)
        flag = 0;
    if (flag == 1)
      return qpow(i, (mod - 1) / k);
  }
}

struct Matrix {
  LL a[3][3];
  Matrix operator * (const Matrix &A) {
    Matrix re;
    for (RI i = 0; i < n; ++i)
      for (RI j = 0; j < n; ++j) {
        re.a[i][j] = 0;
        for (RI k = 0; k < n; ++k)
          re.a[i][j] += a[i][k] * A.a[k][j];
        re.a[i][j] %= mod;
      }
    return re;
  }
  Matrix operator * (const LL &A) {
    Matrix re;
    for (RI i = 0; i < n; ++i)
      for (RI j = 0; j < n; ++j)
        re.a[i][j] = a[i][j] * A % mod;
    return re;
  }
  Matrix operator + (const Matrix &A) {
    Matrix re;
    for (RI i = 0; i < n; ++i)
      for (RI j = 0; j < n; ++j)
        re.a[i][j] = (a[i][j] + A.a[i][j]) % mod;
    return re;
  }
} E, I, F;
Matrix Matrixqpow(Matrix a, int k) {
  Matrix re; re = I;
  for (; k; k >>= 1, a = a * a)
    if (k & 1) re = re * a;
  return re;
}
LL A[MAXN], B[MAXN];

namespace FuckFuckTmd {
  double const PI = acos(-1.0);
  int const M = 32767;
  int r[MAXN];
  struct Complex {
    double a, b;
    Complex (double _a = 0, double _b = 0) { a = _a, b = _b; }
    Complex operator + (const Complex &A) { return Complex(a + A.a, b + A.b); }
    Complex operator - (const Complex &A) { return Complex(a - A.a, b - A.b); }
    Complex operator * (const Complex &A) { return Complex(a * A.a - b * A.b, a * A.b + b * A.a); }
  } omega[MAXN], A[MAXN], B[MAXN], C[MAXN], D[MAXN], F[MAXN], G[MAXN], H[MAXN];
  LL RE[MAXN];
  void FFT(Complex *a, int len, int op) {
    for (RI i = 0; i < len; ++i)
      if (r[i] > i)
        swap(a[i], a[r[i]]);
    for (RI i = 2; i <= len; i <<= 1) {
      int wn = len / i;
      for (RI j = 0; j < len; j += i) {
        int w = 0;
        for (RI k = j; k < j + i / 2; ++k) {
          Complex x = a[k], y = a[k + i / 2] * omega[w];
          a[k] = x + y;
          a[k + i / 2] = x - y;
          w = (w + op * wn + len) % len;
        }
      }
    }
  }
  void FUCKFFT(LL *a, int lena, LL *b, int lenb) {
    for (RI i = 0; i < lena; ++i) {
      A[i].a = (a[i] >> 15) & M;
      B[i].a = a[i] & M;
    }
    for (RI i = 0; i < lenb; ++i) {
      C[i].a = (b[i] >> 15) & M;
      D[i].a = b[i] & M;
    }
    int len = 1, cnt = 0;
    while (len < lena + lenb - 1)
      len <<= 1, ++cnt;
    for (RI i = 0; i < len; ++i)
      r[i] = (r[i >> 1] >> 1) | ((i & 1) << (cnt - 1));
    for (RI i = 0; i < len; ++i) {
      double deg = 2.0 * PI * i / len;
      omega[i] = Complex(cos(deg), sin(deg));
    }
    FFT(A, len, 1), FFT(B, len, 1);
    FFT(C, len, 1), FFT(D, len, 1);
    for (RI i = 0; i < len; ++i) {
      F[i] = A[i] * C[i];
      G[i] = A[i] * D[i] + B[i] * C[i];
      H[i] = B[i] * D[i];
    }
    FFT(F, len, -1), FFT(G, len, -1), FFT(H, len, -1);
    for (RI i = 0; i < len; ++i) {
      RE[i] = LL(F[i].a / len + 0.5) % mod * (M + 1) % mod * (M + 1) % mod;
      RE[i] = (RE[i] + LL(G[i].a / len + 0.5) % mod * (M + 1) % mod) % mod;
      RE[i] = (RE[i] + LL(H[i].a / len + 0.5) % mod) % mod;
      RE[i] = (RE[i] % mod + mod) % mod;
    }
  }
}

//LL ans[MAXN];
int main() {
  
#ifdef LOCAL
  FILEIO("a");
#endif

  cin >> n >> k >> L >> x >> y >> mod; --x, --y;
  for (RI i = 0; i < n; ++i)
    for (RI j = 0; j < n; ++j) {
      cin >> E.a[i][j];
      I.a[i][j] = (i == j);
    }
  LL wn = Getwn(), w = 1;
  int lenF = k, lenG = 2 * k - 1;
  for (RI i = 0; i < lenF; ++i) {
    F = E * w + I;
    F = Matrixqpow(F, L);
    A[i] = F.a[x][y];
    A[i] = A[i] * qpow(wn, 1ll * i * (i - 1) / 2) % mod;
    w = w * wn % mod;
  }
  for (RI i = 0; i < lenG; ++i)
    B[i] = qpow(wn, k - 1ll * i * (i - 1) / 2 % k);
  // for (RI i = 0; i < lenF; ++i)
  //   for (RI j = i; j < lenG; ++j)
  //     ans[j - i] = (ans[j - i] + A[i] * B[j] % mod) % mod;
  for (RI i = 0, j = lenG - 1; i < j; ++i, --j)
    swap(B[i], B[j]);
  FuckFuckTmd :: FUCKFFT(A, lenF, B, lenG);
  LL *ans = FuckFuckTmd :: RE;
  for (RI i = 0; i < k; ++i) {
    LL val = ans[lenG - 1 - i];
    val = val * qpow(k, mod - 2) % mod * qpow(wn, 1ll * i * (i - 1) / 2) % mod;
    cout << val << endl;
  }

  return 0;
}

// created by Daniel yuan
/*
     ________
    /        \
   / /      \ \
  / /        \ \
  \            /
   \  ______  /
    \________/
*/