「HNOI2019」白兔之舞

题目链接

我们先来考虑，如果给定 \(m\)（长度），怎么求方案数。

显然我们可以把横着（二元组的第一维）的单独考虑，即我们要在 \(L\) 个横坐标内选出 \(m\) 个经过的横坐标，这部分的方案为 \(\dbinom{L}{m}\)；然后考虑纵坐标，当 \(n=1\) 的时候显然是 \(w_{1,1}^m\)，当 \(n>1\) 时，如果上一步在 \(a\)，发现会给下一步的 \(1-n\) 分别贡献 \(w_{a,1},w_{a,2},\dots,w_{a,n}\)，发现这是一个矩阵乘法的形式，记 \(W\) 为 \[ \begin{bmatrix} w_{1,1} & \dots & w_{1,n}\\ \vdots & \ddots & \vdots \\ w_{n,1} & \dots & w_{n,n} \end{bmatrix} \] 即第二维的邻接矩阵，为了防止变量重复，我们记题目中的 \(x,y\) 为 \(a,b\)，那么 \(\dbinom{L}{m}W^m_{a,b}\) 就是答案。

然后我们考虑用多项式来表示这个式子，即我们构造一个系数为矩阵的多项式： \[ F(x)=(I+Wx)^L \] 其中 \(I\) 是单位矩阵，则 \(\left([x^m]F(x)\right)_{a,b}\) 就是答案。

那么我们考虑题目要求的是所有模 \(k\) 等于 \(t\) 的 \(m\) 的答案之和，即 \[ \sum\limits_{m\equiv t\pmod k}\left([x^m]F(x)\right)_{a,b} \] 那么我们就是要求 \[ ans_t=\sum\limits_{m\equiv t\pmod k}\left([x^m]F(x)\right)_{a,b} \] 根据单位根反演的简单推广，可以得到 \[ \begin{aligned} ans_t&=\sum\limits_{m\equiv t\pmod k}\left([x^m]F(x)\right)_{a,b}\\ &=\frac{1}{k}\sum\limits_{j=0}^{k-1}w_k^{-tj}F(w_k^j)_{a,b} \end{aligned} \] 推导过程可以看这里，注意单位根反演对矩阵也是成立的（似乎只需要满足数乘对加法有分配律就可以）。

显然 \(F(w_k^j)_{a,b}\) 可以通过矩阵快速幂来求（把邻接矩阵乘上 \(w_k^j\) 再加上单位矩阵的 \(L\) 次方），现在的问题就是我们要对每个 \(0\leqslant t\leqslant k-1\) 都求出 \(ans_t\)。

发现这是一个任意基 IDFT 的形式，并且不一定存在 \(2k\) 次单位根，我们考虑使用 Bluestein's Algorithm 的形式二，即把 \(-tj\) 拆成 \(\dbinom{t}{2}+\dbinom{j}{2}-\dbinom{t+j}{2}\)，那么带入并化简： \[ \begin{aligned} ans_t&=\frac{1}{k}\sum\limits_{j=0}^{k-1}w_k^{-tj}F(w_k^j)_{a,b}\\ &=\frac{1}{k}\sum\limits_{j=0}^{k-1}w_k^{\binom{t}{2}+\binom{j}{2}-\binom{t+j}{2}}F(w_k^j)_{a,b}\\ &=\frac{w_k^{\binom{t}{2}}}{k}\sum\limits_{j=0}^{k-1}w_k^{\binom{j}{2}-\binom{t+j}{2}}F(w_k^j)_{a,b}\\ \end{aligned} \] 设 \(A(x)=\sum\limits_iw_k^{\binom{i}{2}}F(w_k^j)_{a,b}x^i\)，\(B(x)=\sum\limits_iw_k^{-\binom{i}{2}}x^i\)后面是一个差一定的卷积形式，可以通过反转 \(B\) 求解。

由于这题不一定是 NTT 模数，所以需要自己找到原根然后写一个 MTT。

代码：

#include<bits/stdc++.h>
using namespace std;
int F[400005],G[400005],H[400005],g,n,k,l,a,b,p;
namespace poly
{
    long double const pi=acos(-1);
    struct comp
    {
        long double r,i;
        comp(){r=i=0;}
        comp(long double x,long double y){r=x,i=y;}
        comp conj(){return comp(r,-i);}
        friend comp operator +(comp x,comp y){return comp(x.r+y.r,x.i+y.i);}
        friend comp operator -(comp x,comp y){return comp(x.r-y.r,x.i-y.i);}
        friend comp operator *(comp x,comp y){return comp(x.r*y.r-x.i*y.i,x.i*y.r+x.r*y.i);}
    };
    typedef long long ll;
    int r[400005];
    comp a[400005],b[400005],c[400005],d[400005];
    void fft(comp *f,int n,int op)
    {
        for(int i=1;i<n;i++)r[i]=(r[i>>1]>>1)+((i&1)?(n>>1):0);
        for(int i=1;i<n;i++)if(i<r[i])swap(f[i],f[r[i]]);
        for(int len=2;len<=n;len<<=1)
        {
            int q=len>>1;
            comp wn=comp(cos(pi/q),op*sin(pi/q));
            for(int i=0;i<n;i+=len)
            {
                comp w=comp(1,0);
                for(int j=i;j<i+q;j++,w=w*wn)
                {
                    comp d=f[j+q]*w;
                    f[j+q]=f[j]-d;
                    f[j]=f[j]+d;
                }
            }
        }
    }
    void mtt(int *f,int *g,int *h,int n,int p)
    {
        for(int i=0;i<n;i++)
            a[i].r=(f[i]>>15),a[i].i=(f[i]&32767),
            c[i].r=(g[i]>>15),c[i].i=(g[i]&32767);
        fft(a,n,1),fft(c,n,1);
        for(int i=1;i<n;i++)b[i]=a[n-i].conj();
        b[0]=a[0].conj();
        for(int i=1;i<n;i++)d[i]=c[n-i].conj();
        d[0]=c[0].conj();
        for(int i=0;i<n;i++)
        {
            comp 
            aa=(a[i]+b[i])*comp(0.5,0),
            bb=(a[i]-b[i])*comp(0,-0.5),
            cc=(c[i]+d[i])*comp(0.5,0),
            dd=(c[i]-d[i])*comp(0,-0.5);
            a[i]=aa*cc+comp(0,1)*(aa*dd+bb*cc),b[i]=bb*dd;
        }
        fft(a,n,-1),fft(b,n,-1);
        for(int i=0;i<n;i++)
        {
            int 
            aa=(ll)(a[i].r/n+0.5)%p,
            bb=(ll)(a[i].i/n+0.5)%p,
            cc=(ll)(b[i].r/n+0.5)%p;
            h[i]=((1ll*aa*(1<<30)+1ll*bb*(1<<15)+cc)%p+p)%p;
        }
    }
}
using poly::mtt;
int pw(int x,int y)
{
    int res=1;
    while(y)
    {
        if(y&1)res=1ll*res*x%p;
        x=1ll*x*x%p;
        y>>=1;
    }
    return res;
}
inline int mod(int x){return x>=p?x-p:x;}
int get(int x)
{
    vector<int>v;
    int t=p-1;
    for(int i=2;i*i<=t;i++)
    {
        if(!(t%i))
        {
            v.push_back(i);
            while(!(t%i))t/=i;
        }
    }
    for(int i=2;i<=p-1;i++)
    {
        bool flag=1;
        for(int j=0;j<v.size();j++)
            if(pw(i,(p-1)/v[j])==1){flag=0;break;}
        if(flag)return i;
    }
}
struct matrix
{
    int s[4][4];
    matrix(){memset(s,0,sizeof(s));}
    friend matrix operator *(matrix a,matrix b)
    {
        matrix tmp;
        for(int i=1;i<=n;i++)
            for(int j=1;j<=n;j++)
                for(int k=1;k<=n;k++)
                    tmp.s[i][j]=mod(tmp.s[i][j]+1ll*a.s[i][k]*b.s[k][j]%p);
        return tmp;
    }
    friend matrix operator *(int a,matrix b)
    {
        for(int i=1;i<=n;i++)
            for(int j=1;j<=n;j++)
                b.s[i][j]=1ll*b.s[i][j]*a%p;
        return b;
    }
    friend matrix operator +(matrix a,matrix b)
    {
        for(int i=1;i<=n;i++)
            for(int j=1;j<=n;j++)
                b.s[i][j]=mod(a.s[i][j]+b.s[i][j]);
        return b;
    }
}W,I;
matrix pw2(matrix x,int y)
{
    matrix res=I;
    while(y)
    {
        if(y&1)res=res*x;
        x=x*x;
        y>>=1;
    }
    return res;
}
int main()
{
    scanf("%d%d%d%d%d%d",&n,&k,&l,&a,&b,&p);
    g=get(p);
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            scanf("%d",&W.s[i][j]);
    for(int i=1;i<=n;i++)I.s[i][i]=1;
    int wk=pw(g,(p-1)/k),w=1;
    for(int i=0;i<k;i++,w=1ll*w*wk%p)
        F[i]=1ll*pw(wk,(1ll*i*(i-1)/2)%k)*pw2(w*W+I,l).s[a][b]%p;
    for(int i=0;i<(k<<1);i++)
        G[i]=1ll*pw(wk,(p-1-(1ll*i*(i-1)/2%k)));
    reverse(G,G+(k<<1));
    int lim=1,inv=pw(k,p-2);
    while(lim<(k<<1)+1)lim<<=1;
    mtt(F,G,H,lim,p);
    reverse(H,H+(k<<1));
    for(int i=0;i<k;i++)printf("%lld\n",1ll*inv*H[i]%p*pw(wk,(1ll*i*(i-1)/2)%k)%p);
    return 0;
}