[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+;

struct mat{ll a[][];}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=;

    while(b)

    {

        if(b&)ret=ret*a%mod;

        a=a*a%mod,b>>=;

    }

    return ret;

}

mat operator*(mat a,mat b)

{

    mat c;

    for(int i=;i<=n;i++)

    for(int j=;j<=n;j++)

    {

        c.a[i][j]=;

        for(int k=;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={};

    for(int i=;i<=n;i++)

    for(int j=;j<=n;j++)

    if(i==j)ret.a[i][j]=;

    while(b)

    {

        if(b&)ret=ret*a;

        a=a*a,b>>=;

    }

    return ret;

}

ll getrt(ll x)

{

    vector<int>vec;vec.clear();

    for(ll i=;i*i<=(x-);i++)if((x-)%i==)vec.push_back(i),vec.push_back((x-)/i);

    for(ll i=;;i++)

    {

        bool c=;

        for(unsigned j=;j<vec.size()&&c;j++)if(qpow(i,vec[j])==)c=;

        if(c)return i;

    }

}

ll c1(ll x){return x*(x-)/;}

void FFT(cp*a,int f)

{

    for(int i=;i<nn;i++)if(i<R[i])swap(a[i],a[R[i]]);

    for(int i=;i<nn;i*=)

    {

        cp wn(cosl(pi/i),sinl(pi/i));

        for(int j=;j<nn;j+=i*)

        {

            cp w(,),x,y;

            for(int k=;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==)return;

    reverse(a+,a+nn);

    for(int i=;i<nn;i++)a[i]/=nn;

}

int main()

{

    cin>>n>>k>>l>>X>>Y>>mod;

    for(int i=;i<=n;i++)for(int j=;j<=n;j++)cin>>A1.a[i][j];

    w=qpow(getrt(mod),(mod-)/k);

    for(int j=;j<k;j++)

    {

        for(int i=;i<=n;i++)

        for(int p=;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=;j<k;j++)a2[j]=a2[j]*qpow(w,c1(j))%mod;

    for(int j=;j<=k*;j++)a1[j]=qpow(w,c1(j)*(mod-));

    reverse(a1,a1+k*+);

    int m=k*,L=;

    for(int i=;i<m;i++)A[i]=a2[i]/,B[i]=a2[i]%;

    for(int i=;i<m;i++)C[i]=a1[i]/,D[i]=a1[i]%;

    for(nn=;nn<=m;nn<<=)L++;

    for(int i=;i<nn;i++)R[i]=(R[i>>]>>)|((i&)<<L-);

    FFT(A,),FFT(B,),FFT(C,),FFT(D,);

    for(int i=;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,-),FFT(F,-),FFT(G,-);

    for(int i=;i<nn;i++)

    ans[i]=((ll)(G[i].real()+0.5)+(((ll)(F[i].real()+0.5)%mod)<<)+((ll)(E[i].real()+0.5)%mod<<))%mod;

    for(ll i=;i<k;i++)

    {

        ans[k*-i]=(ans[k*-i]*qpow(k,mod-)%mod)*qpow(w,c1(i))%mod;

        printf("%lld\n",ans[k*-i]);

    }

}