loj #3058. 「HNOI2019」白兔之舞

题意：

自己看

题解：

观察一下，发现行和列的方案是无关的，只要乘起来就好。所以走$m$步结束的方案数就是：

$W$表示第二维的转移矩阵（输入的那个）。现在要求$m \mod k = t$的答案，考虑单位根反演：

这个$\omega^{-it}$不是很好搞。考虑把$it$换成$\binom{i + t}{2} - \binom{i}{2} - \binom{t}{2}$：

这就是一个卷积的形式了，用任意模数fft即可。

代码：

#include <bits/stdc++.h>
using namespace std;

const double pi = acos(-1);
int n, K, l, x, y, mod, qmod, N, a[262144], b[262144], c[262144];
struct matrix
{
    int x, y;
    long long a[3][3];
    matrix() { memset(a, 0, sizeof(a)); }
    long long *operator[](int x) { return a[x]; }
} I, W;
struct Complex
{
    double real, imag;
    Complex(double real = 0, double imag = 0) : real(real), imag(imag) {}
} w[262144], a1[262144], a2[262144], b1[262144], b2[262144], c1[262144], c2[262144], c3[262144];
Complex operator+(Complex x, Complex y) { return Complex(x.real + y.real, x.imag + y.imag); }
Complex operator-(Complex x, Complex y) { return Complex(x.real - y.real, x.imag - y.imag); }
Complex operator*(Complex x, Complex y) { return Complex(x.real * y.real - x.imag * y.imag, x.real * y.imag + x.imag * y.real); }
int g[65536];
char buf[1 << 20], *now = buf;

void flush()
{
    fwrite(buf, 1, now - buf, stdout);
    now = buf;
}
void write(int x)
{
    static int a[10];
    int num = 0;
    if (x == 0) a[num++] = 0;
    while (x)
    {
        a[num++] = x % 10;
        x /= 10;
    }
    for (int i = num - 1; i >= 0; i--)
    {
        *(now++) = a[i] + '0';
        if (now == buf + (1 << 20)) flush();
    }
    *(now++) = '\n';
    if (now == buf + (1 << 20)) flush();
}
matrix operator*(matrix x, matrix y)
{
    matrix ans;
    ans.x = x.x;
    ans.y = y.y;
    for (int i = 0; i < x.x; i++)
        for (int j = 0; j < y.y; j++)
            for (int k = 0; k < y.x; k++)
                (ans[i][j] += x[i][k] * y[k][j] % mod) %= mod;
    return ans;
}
matrix operator+(matrix x, matrix y)
{
    for (int i = 0; i < x.x; i++)
        for (int j = 0; j < x.y; j++)
            (x[i][j] += y[i][j]) %= mod;
    return x;
}
matrix operator*(matrix x, long long y)
{
    for (int i = 0; i < x.x; i++)
        for (int j = 0; j < x.y; j++)
            (x[i][j] *= y) %= mod;
    return x;
}
matrix Pow(matrix x, int y)
{
    matrix ans = I;
    for (; y; y >>= 1, x = x * x)
        if (y & 1) ans = ans * x;
    return ans;
}
void bit_reverse(Complex *x, int n)
{
    for (int i = 0, j = 0; i < n; i++)
    {
        if (i > j) swap(x[i], x[j]);
        for (int k = n >> 1; (j ^= k) < k; k >>= 1);
    }
}
void fft(Complex *x, int n, int op)
{
    bit_reverse(x, n);
    for (int i = 0; i < n; i++) w[i].imag *= op;
    for (int i = 2; i <= n; i <<= 1)
    {
        int m = i >> 1;
        for (int j = 0; j < n; j += i)
            for (int k = 0; k < m; k++)
            {
                Complex hh = x[j + m + k] * w[n / i * k];
                x[j + m + k] = x[j + k] - hh;
                x[j + k] = x[j + k] + hh;
            }
    }
    if (op == -1)
        for (int i = 0; i < n; i++)
        {
            w[i].imag *= -1;
            x[i].real = x[i].real / n;
        }
}
int Pow(int x, int y)
{
    int ans = 1;
    for (; y; y >>= 1, x = 1LL * x * x % mod)
        if (y & 1) ans = 1LL * ans * x % mod;
    return ans;
}
int get_g()
{
    static vector<int> p;
    int hh = mod - 1;
    for (int i = 2; i * i <= hh; i++)
        if (hh % i == 0)
        {
            p.push_back(i);
            while (hh % i == 0) hh /= i;
        }
    if (hh > 1) p.push_back(hh);
    int g = 2;
    while (1)
    {
        bool tf = false;
        for (int j = 0; j < p.size(); j++)
            if (Pow(g, (mod - 1) / p[j]) == 1)
            {
                tf = true;
                break;
            }
        if (!tf) return g;
        g++;
    }
}
int main()
{
    scanf("%d%d%d%d%d%d", &n, &K, &l, &x, &y, &mod);
    I.x = I.y = W.x = W.y = n;
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            scanf("%lld", &W[i][j]);
    for (int i = 0; i < n; i++) I[i][i] = 1;
    int gk = Pow(get_g(), (mod - 1) / K);
    g[0] = 1;
    for (int i = 1; i < K; i++) g[i] = 1LL * g[i - 1] * gk % mod;
    for (int i = 0; i < K; i++)
    {
        a[i] = 1LL * g[1LL * i * (i - 1) / 2 % K] * Pow(W * g[i] + I, l).a[x - 1][y - 1] % mod;
    }
    for (int i = 0; i < (K << 1); i++)
        b[K + K - i - 1] = g[(K - (1LL * i * (i - 1) / 2 % K)) % K];
    qmod = sqrt(mod);
    N = 1;
    while (N < (K << 2)) N <<= 1;
    for (int i = 0; i < N; i++)
    {
        w[i].real = cos(2 * pi / N * i);
        w[i].imag = sin(2 * pi / N * i);
    }
    for (int i = 0; i < (K << 1); i++)
    {
        a1[i].real = a[i] / qmod;
        a2[i].real = a[i] % qmod;
        b1[i].real = b[i] / qmod;
        b2[i].real = b[i] % qmod;
    }
    fft(a1, N, 1);
    fft(a2, N, 1);
    fft(b1, N, 1);
    fft(b2, N, 1);
    for (int i = 0; i < N; i++)
    {
        c1[i] = a1[i] * b1[i];
        c2[i] = a1[i] * b2[i] + a2[i] * b1[i];
        c3[i] = a2[i] * b2[i];
    }
    fft(c1, N, -1);
    fft(c2, N, -1);
    fft(c3, N, -1);
    for (int i = 0; i < N; i++)
        c[i] = (llround(c1[i].real) * qmod % mod * qmod % mod + llround(c2[i].real) * qmod % mod + llround(c3[i].real) % mod) % mod;
    int invK = Pow(K, mod - 2);
    for (int i = (K << 1) - 1; i >= K; i--)
        // printf("%lld\n", 1LL * c[i] * invK % mod * g[(1LL * (K + K - i - 1) * (K + K - i - 2) / 2 % K) % K] % mod);
        write(1LL * c[i] * invK % mod * g[(1LL * (K + K - i - 1) * (K + K - i - 2) / 2 % K) % K] % mod);
    flush();
}