[HNOI2019]白兔之舞（矩阵快速幂+单位根反演）

非常抱歉，这篇文章鸽了。真的没时间写了。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
#define cp complex<long double>
#define pi acosl(-1)
const int N=6e5+7;
struct mat{ll a[5][5];}A1,A2;
ll n,k,l,X,Y,mod,w,nn,ans[N],a1[N],a2[N],R[N];
cp A[N],B[N],C[N],D[N],E[N],F[N],G[N];
ll qpow(ll a,ll b)
{
    ll ret=1;
    while(b)
    {
        if(b&1)ret=ret*a%mod;
        a=a*a%mod,b>>=1;
    }
    return ret;
}
mat operator*(mat a,mat b)
{
    mat c;
    for(int i=1;i<=n;i++)
    for(int j=1;j<=n;j++)
    {
        c.a[i][j]=0;
        for(int k=1;k<=n;k++)c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j])%mod;
    }
    return c;
}
mat qpow(mat a,ll b)
{
    mat ret={0};
    for(int i=1;i<=n;i++)
    for(int j=1;j<=n;j++)
    if(i==j)ret.a[i][j]=1;
    while(b)
    {
        if(b&1)ret=ret*a;
        a=a*a,b>>=1;
    }
    return ret;
}
ll getrt(ll x)
{
    vector<int>vec;vec.clear();
    for(ll i=2;i*i<=(x-1);i++)if((x-1)%i==0)vec.push_back(i),vec.push_back((x-1)/i);
    for(ll i=2;;i++)
    {
        bool c=1;
        for(unsigned j=0;j<vec.size()&&c;j++)if(qpow(i,vec[j])==1)c=0;
        if(c)return i;
    }
}
ll c1(ll x){return x*(x-1)/2;}
void FFT(cp*a,int f)
{
    for(int i=0;i<nn;i++)if(i<R[i])swap(a[i],a[R[i]]);
    for(int i=1;i<nn;i*=2)
    {
        cp wn(cosl(pi/i),sinl(pi/i));
        for(int j=0;j<nn;j+=i*2)
        {
            cp w(1,0),x,y;
            for(int k=0;k<i;k++,w=w*wn)
            x=a[j+k],y=a[j+k+i]*w,a[j+k]=x+y,a[j+k+i]=x-y;
        }
    }
    if(f==1)return;
    reverse(a+1,a+nn);
    for(int i=0;i<nn;i++)a[i]/=nn;
}
int main()
{
    cin>>n>>k>>l>>X>>Y>>mod;
    for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)cin>>A1.a[i][j];
    w=qpow(getrt(mod),(mod-1)/k);
    for(int j=0;j<k;j++)
    {
        for(int i=1;i<=n;i++)
        for(int p=1;p<=n;p++)
        {
            A2.a[i][p]=A1.a[i][p]*qpow(w,j)%mod;
            if(i==p)A2.a[i][p]++;
        }
        A2=qpow(A2,l),a2[j]=A2.a[X][Y];
    }
    for(int j=0;j<k;j++)a2[j]=a2[j]*qpow(w,c1(j))%mod;
    for(int j=0;j<=k*2;j++)a1[j]=qpow(w,c1(j)*(mod-2));
    reverse(a1,a1+k*2+1);
    int m=k*3,L=0;
    for(int i=0;i<m;i++)A[i]=a2[i]/32768,B[i]=a2[i]%32768;
    for(int i=0;i<m;i++)C[i]=a1[i]/32768,D[i]=a1[i]%32768;
    for(nn=1;nn<=m;nn<<=1)L++;
    for(int i=0;i<nn;i++)R[i]=(R[i>>1]>>1)|((i&1)<<L-1);
    FFT(A,1),FFT(B,1),FFT(C,1),FFT(D,1);
    for(int i=0;i<nn;i++)E[i]=A[i]*C[i],F[i]=A[i]*D[i]+B[i]*C[i],G[i]=B[i]*D[i];
    FFT(E,-1),FFT(F,-1),FFT(G,-1);
    for(int i=0;i<nn;i++)
    ans[i]=((ll)(G[i].real()+0.5)+(((ll)(F[i].real()+0.5)%mod)<<15)+((ll)(E[i].real()+0.5)%mod<<30))%mod;
    for(ll i=0;i<k;i++)
    {
        ans[k*2-i]=(ans[k*2-i]*qpow(k,mod-2)%mod)*qpow(w,c1(i))%mod;
        printf("%lld
",ans[k*2-i]);
    }
}