[HNOI 2019] 白兔之舞

LOJ #3058.

看着这个题感觉就有点不可做……所以我们先看$n = 1$的情况

设$c_i$为走$i$步,从$(0, x)$走到一个坐标第二维为$y$的顶点的方案数,那么显然有
$$c_i = {n\choose i}\cdot {w_{1, 1}}^i$$

考虑枚举走的步数统计答案,答案即为
$$\begin{split} ans_t
&= \sum \limits_{i=0}^Lc_i\cdot [k | i-t]\\
&= \sum \limits_{i=0}^L{n\choose i}\cdot {w_{1, 1}}^i\cdot \frac{1}{k}\cdot \sum \limits_{j=0}^{k-1}{\omega_k}^{j(i-t)}\\
&= \frac{1}{k}\cdot \sum \limits_{j=0}^{k-1}{\omega_{k}}^{-jt}\cdot \sum \limits_{i=0}^L{n\choose i}\cdot {w_{1, 1}}^i\cdot {\omega_k}^{ji}\\
&=\frac{1}{k}\cdot \sum \limits_{j=0}^{k-1}{\omega_k}^{-jt}\cdot({w_{1, 1}}\cdot {\omega_k}^{j} + 1)^L\\
&=\frac{1}{k}\cdot \sum \limits_{j=0}^{k-1}w_k^{-{j+t\choose 2}+{j\choose 2}+{t\choose 2}}\cdot ({w_{1, 1}}\cdot {\omega_k}^{j} + 1)^L\\
&=\frac{{\omega_{k}}^{t\choose 2}}{k}\cdot \sum \limits_{j =0}^{k-1}{\omega_k}^{-{j+t\choose 2}}\cdot {\omega_k}^{j\choose 2}\cdot (w_{1, 1}\cdot{\omega_k}^j+1)^L
\end{split}$$

其中第一行到第二行套了一个单位根反演,第四行到第五行运用了一个$ij={i+j\choose 2}-{i\choose 2}-{j\choose 2}$的高端操作,然后这个最后一步显然可以卷积了,设$f_{k-i-1}={\omega_k}^{-{i\choose 2}},g_{i} = {\omega_k}^{i\choose 2}\cdot({\omega_k}^j+1)^L$即可

至于$n\not = 1$的情况,将$w_{1, 1}$改为输入的$w$矩阵,$1$改为单位矩阵$I$就可以了

#include <iostream>
#include <cstdio>
#include <vector>
#include <cmath>
#include <cstring>

namespace IO{
#define MX (1 << 21)
    char ibuf[MX], obuf[MX], *iS, *iT, *oS = obuf, *oT = obuf + MX, c;
    inline char getc(){ return getchar();
        if(iS == iT) iT = (iS = ibuf) + fread(ibuf, 1, MX, stdin);
        return *iS++;
    }
    inline void flush(){
        fwrite(obuf, 1, oS - obuf, stdout);
        oS = obuf;
    }
    inline void putc(char _){
        if(oS == oT) flush(); *oS++ = _;
    }
    struct autof{ ~autof(){flush();} }ff;
    template <typename T> inline void in(T &_){
        c = getc(), _ = 0; int fl = 1;
        while(c < '0' || c > '9') 
            fl = c == '-' ? -1 : fl, c = getc();
        while(c >= '0' && c <= '9') 
            _ = _ * 10 + c - '0', c = getc(); 
        _ *= fl;
    }
    template <typename T> void out(T _){
        if(_ < 0) putc('-'), out(-_);
        else{
            if(_ / 10) out(_ / 10); 
            putc('0' + _ % 10);
        }
    }
#undef MX
}
using IO::in;
using IO::out;
using IO::putc;

#define pb push_back
#define MAX_K 65546
typedef long long ll;

int N, K, L, X, Y, P;

int ksm(int x, int y){
    int z = 1;
    for(; y; y >>= 1, x = 1ll * x * x % P)
        if(y & 1) z = 1ll * x * z % P;
    return z;
}

namespace ROOT{
    int w[MAX_K];

    int o(int x){
        x = 1ll * x * (x-1) / 2 % K;
        return w[x < 0 ? x + K : x];
    }

    std::vector<int> get_pri(int x){
        std::vector<int> res;
        for(int i = 2; i * i <= x; i++) if(x % i == 0){
            while(x % i == 0) x /= i;
            res.pb(i);
        }
        if(x > 1) res.pb(x);
        return res;
    }

    int get_rt(int x){
        std::vector<int> pri = get_pri(x - 1);
        for(int i = 2; i < x; i++){
            for(auto j: pri) if(ksm(i, (x - 1) / j) == 1)
                goto qwq;
            return i;
            qwq:;
        }
        return -1;
    }

    void pre(){
        w[0] = 1, w[1] = ksm(get_rt(P), (P - 1) / K);
        for(int i = 2; i < K; i++)
            w[i] = 1ll * w[i-1] * w[1] % P;
    }
}
using ROOT::o;

struct Matrix{
    int a[4][4];
    Matrix(){
        memset(a, 0, sizeof(a));
    }
    int *operator [](int x){
        return a[x];
    }
}A, I;

Matrix operator +(Matrix x, Matrix y){
    for(int i = 1; i <= N; i++)
        for(int j = 1; j <= N; j++)
            x[i][j] = (x[i][j] + y[i][j]) % P;
    return x;
}

Matrix operator *(Matrix x, Matrix y){
    Matrix z;
    for(int i = 1; i <= N; i++)
        for(int j = 1; j <= N; j++)
            for(int k = 1; k <= N; k++)
                z[i][j] = (z[i][j] + 1ll * x[i][k] * y[k][j] % P) % P;
    return z;
}

Matrix operator *(int x, Matrix y){
    for(int i = 1; i <= N; i++)
        for(int j = 1; j <= N; j++)
            y[i][j] = 1ll * x * y[i][j] % P;
    return y;
}

Matrix ksm(Matrix x, int y){
    Matrix z;
    for(int i = 1; i <= N; i++) z[i][i] = 1;
    for(; y; y >>= 1, x = x * x)
        if(y & 1) z = z * x;
    return z;
}

struct Poly : std::vector<int>{};

namespace MTT{
#define MX 400010
    const long double Pi = std::acos(-1);
    int rev[MX];
    struct Complex{
        long double x, y;
        Complex(){}
        Complex(long double x, long double y) : x(x), y(y){}
    }a[MX], b[MX], c[MX], d[MX], e[MX], f[MX], g[MX], h[MX], w[MX], iw[MX];
    Complex operator +(Complex x, Complex y){
        return Complex(x.x + y.x, x.y + y.y);
    }
    Complex operator -(Complex x, Complex y){
        return Complex(x.x - y.x, x.y - y.y);
    }
    Complex operator *(Complex x, Complex y){
        return Complex(x.x * y.x - x.y * y.y, x.x * y.y + x.y * y.x);
    }

    void fft(Complex *P, int l, int opt){
        for(int i = 0; i < l; i++) if(rev[i] > i) 
        std::swap(P[rev[i]], P[i]);
        for(int i = 1; i < l; i <<= 1){
            for(int j = 0, p = i << 1; j < l; j += p)
                for(int k = 0; k < i; k++){
                    Complex x = P[j + k], y = P[i + j + k] * (opt > 0 ? w[i + k] : iw[i + k]);
                    P[j + k] = x + y, P[i + j + k] = x - y;
                }
        }
        if(opt < 0)
            for(int i = 0; i < l; i++) P[i].x /= l;
    }

    Poly mul(Poly x, Poly y){
        int l = 1;
        while(l < x.size() + y.size() - 1) l <<= 1;
        l >>= 1;
        for(int i = 0; i < l; i++) 
            w[i + l] = Complex(std::cos(Pi * i / l), std::sin(Pi * i / l));
        for(int i = 0; i < l; i++) 
            iw[i + l] = Complex(std::cos(Pi * i / l), -std::sin(Pi * i / l));
        for(int i = l - 1; ~i; i--) 
            w[i] = w[i << 1], iw[i] = iw[i << 1];
        l <<= 1;
        for(int i = 0; i < l; i++)
            rev[i] = (rev[i>>1] >> 1) | ((i&1) ? l>>1 : 0);
        for(int i = 0; i < x.size(); i++){
            a[i] = Complex(x[i] >> 15, 0);
            b[i] = Complex(x[i] & 32767, 0);
        }
        for(int i = 0; i < y.size(); i++){
            c[i] = Complex(y[i] >> 15, 0);
            d[i] = Complex(y[i] & 32767, 0);
        }
        for(int i = x.size(); i < l; i++)
            a[i] = Complex(0, 0), b[i] = Complex(0, 0);
        for(int i = y.size(); i < l; i++)
            c[i] = Complex(0, 0), d[i] = Complex(0, 0);
        fft(a, l, 1); fft(b, l, 1); fft(c, l, 1); fft(d, l, 1);
        for(int i = 0; i < l; i++){
            e[i] = a[i] * c[i], f[i] = a[i] * d[i];
            g[i] = b[i] * c[i], h[i] = b[i] * d[i];
        } 
        fft(e, l, -1); fft(f, l, -1); fft(g, l, -1); fft(h, l, -1);
        int lim = x.size() + y.size() - 1;
        x.resize(lim);
        ll p30 = (1 << 30) % P, p15 = (1 << 15) % P;
        for(int i = 0; i < lim; i++){
            x[i] = (ll)(round(e[i].x)) % P * p30 % P;
            x[i] = ((ll)(round(f[i].x)) % P * p15 % P+ x[i])%P;
            x[i] = ((ll)(round(g[i].x)) % P * p15 % P+ x[i])%P;
            x[i] = (x[i] + (ll)(round(h[i].x)) % P)%P;
         }
        return x;
    }
#undef MX
}

int main(){
    in(N), in(K), in(L), in(X), in(Y), in(P);
    ROOT::pre();
    for(int i = 1; i <= N; i++)
        for(int j = 1; j <= N; j++)
            in(A[i][j]);
    for(int i = 1; i <= N; i++) I[i][i] = 1;
    Poly f, g;
    g.resize(K); f.resize((K << 1) - 1);
    for(int i = 0; i < K; i++){
        int c = ksm(ROOT::w[i] * A + I, L)[X][Y];
        g[K - i - 1] = 1ll * c * o(i) % P;
    }
    for(int i = 0; i < (K << 1) - 1; i++)
        f[i] = ksm(o(i), P - 2);
    f = MTT::mul(f, g);
    for(int i = 0; i < K; i++)
        printf("%lld\n", 1ll * f[i + K - 1] * ksm(K, P - 2) % P * o(i) % P);
    return 0;
}

发表评论

电子邮件地址不会被公开。 必填项已用*标注