[HNOI2019]白兔之舞(MTT+单位根反演)

[HNOI2019]白兔之舞

\[注:W是一个矩阵,表示题中w[i][j],下列式子中的W做加减法时表示W_{xy}\\ ans_t=\Sigma_{i\ mod\ k=t}^L{L \choose i}W^i\\ =\Sigma_{i=0}^L[k|(i-t)]{L \choose i}W^i\\ =\Sigma_{i=0}^L\frac1k\Sigma_{j=0}^{k-1}(\omega_k^{i-t})^j{L \choose i}W^i\\ =\Sigma_i^L\frac1k\Sigma_{j=0}^{k-1}\omega_k^{ij-tj}{L \choose i}W^i\\ =\frac1k\Sigma_{j=0}^{k-1}\omega_k^{{t \choose 2}+{j \choose 2}-{t+j \choose 2}}\Sigma_{i=0}^L\omega_k^{ij}{L \choose i}W^i\\ 设C(j)=\Sigma_{i=0}^L\omega_k^{ij}{L \choose i}W^i\\ 把它变成二项式的形式:\\ A(j)=\Sigma_{i=0}^L{L \choose i}(W\omega_k^j)^i*1^{L-i}\\ \therefore A(j)=(W\omega_k^j+1)^L\\ \therefore ans_t=\frac{\omega_k^{t \choose 2}}k\Sigma_{j=0}^{k-1}\omega_k^{j \choose 2}A(j)\omega_k^{-{t+j \choose 2}}\\ 这好像是还不是卷积,但是我们可以把它变成卷积\\ 令f(x)=\omega_k^{x \choose 2}\\ 令g(x)=A(x)\omega_k^{-{t+x \choose 2}}\\ 那么按照常规操作,将f倍增,再将g翻转.\\ \therefore ans_t=\frac{\omega_k^{t \choose 2}}k(f*g)(k+t)\\ 然后卷积用MTT算,就没了. \]

#pragma GCC optimize(3)
#include<bits/stdc++.h>
#define N (65536*4)
using namespace std;
const double pi=acos(-1);
int n,k,l,x,y,p;
struct Matrix{
    int a[3][3];
    void format(){
        memset(a,0,sizeof(a));
    }
    friend Matrix operator + (const Matrix &a,const Matrix &b){
        Matrix re;
        for(int i=0;i<3;i++)
        for(int j=0;j<3;j++){
            re.a[i][j]=(a.a[i][j]+b.a[i][j])%p;
        }
        return re;
    }
    friend Matrix operator * (const Matrix &a,const Matrix &b){
        Matrix re;
        re.format();
        for(int i=0;i<3;i++)
        for(int j=0;j<3;j++)
        for(int k=0;k<3;k++){
            re.a[i][j]=(re.a[i][j]+1ll*a.a[i][k]*b.a[k][j])%p;
        }
        return re;
    }
    friend Matrix operator * (const int &a,const Matrix &b){
        Matrix re;
        for(int i=0;i<3;i++)
        for(int j=0;j<3;j++){
            re.a[i][j]=(1ll*a*b.a[i][j])%p;
        }
        return re;
    }
    friend Matrix operator * (const Matrix &a,const int &b){return b*a;}
}I,W;
Matrix Pow(Matrix x,int y){
    Matrix re=I;
    while(y){
        if(y&1)re=re*x;
        y>>=1;
        x=x*x;
    }
    return re;
}
struct Complex{double x,y;};
Complex operator + (Complex a,Complex b){return (Complex){a.x+b.x,a.y+b.y};}
Complex operator - (Complex a,Complex b){return (Complex){a.x-b.x,a.y-b.y};}
Complex operator * (Complex a,Complex b){return (Complex){a.x*b.x-a.y*b.y,a.y*b.x+a.x*b.y};}
Complex a[N],b[N],c[N];

int Pow(int x,int y){
    int re=1;
    while(y){
        if(y&1)re=(1ll*re*x)%p;
        y>>=1;
        x=(1ll*x*x)%p;
    }
    return re;
}

int rev[N];
int pre(int n,int m){
    int high=0,cnt=1;;
    while(cnt<n+m)cnt<<=1,high++;
    for(int i=0;i<cnt;i++)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(high-1));
    return cnt;
}
int root(int num){
    int div[128],sum=0;
    int phi=num-1;
    for(int i=2;i<=phi;i++){
        if(phi%i==0){
            div[++sum]=i;
            while(phi%i==0)phi/=i;
        }
    }
    if(phi!=1)div[++sum]=phi;
    int re=2;
    while(1){
        int flag=1;
        for(int i=1;i<=sum;i++){
            int pw=Pow(re,(num-1)/div[i]);
            if(pw==1){flag=0;break;}
        }
        if(flag)return re;
        re++;
    }
}
Complex w[N];
void fft(int lim,Complex *buf,int dft=1){
    for(int i=0;i<lim;i++)if(i<rev[i])swap(buf[i],buf[rev[i]]);
    for(int len=2;len<=lim;len<<=1){
        for(int s=0;s<lim;s+=len){
            for(int k=s;k<s+len/2;k++){
                Complex x,y;
                x=buf[k],y=w[lim/len*(k-s)]*buf[k+len/2];
                buf[k]=x+y;
                buf[k+len/2]=x-y;
            }
        }
    }
    if(dft==-1)for(int i=0;i<lim;i++)buf[i].x/=lim;
}
Complex A1[N],B1[N],A2[N],B2[N];
void Mul(int lim,int *A,int *B){
    for(int i=0;i<lim;i++)A1[i].x=A[i]>>15,B1[i].x=A[i]&32767;
    for(int i=0;i<lim;i++)A2[i].x=B[i]>>15,B2[i].x=B[i]&32767;
    for(int i=0;i<N;++i)w[i]=(Complex){cos(i*pi*2/lim),sin(i*pi*2/lim)};
    fft(lim,A1);
    fft(lim,B1);
    fft(lim,A2);
    fft(lim,B2);
    for(int i=0;i<lim;i++){
        Complex A1_=A1[i],A2_=A2[i];
        Complex B1_=B1[i],B2_=B2[i];
        A1[i]=A1_*A2_;
        A2[i]=B1_*B2_;
        B1[i]=A1_*B2_;
        B2[i]=A2_*B1_;
    }
    for(int i=0;i<lim;++i)w[i]=(Complex){cos(i*pi*2/lim),-sin(i*pi*2/lim)};
    fft(lim,A1,-1);
    fft(lim,B1,-1);
    fft(lim,A2,-1);
    fft(lim,B2,-1);
    int m=32768;
    for(int i=0;i<lim;i++){
        double A1_=A1[i].x,A2_=A2[i].x;
        double B1_=B1[i].x,B2_=B2[i].x;
        A[i]=((long long)(A1_+0.5)*m%p*m%p+(long long)(B1_+B2_+0.5)*m%p+(long long)(A2_+0.5))%p;
    }
}
int gen[N],f[N],g[N];
int omg(long long x){
    x=((x%k)+k)%k;
    return gen[x];
}
signed main(){
#ifdef KOG_juruo
    freopen("data.in","r",stdin);
    freopen("my.out","w",stdout);
#endif
    scanf("%d%d%d%d%d%d",&n,&k,&l,&x,&y,&p);
    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            scanf("%d",&W.a[i][j]);
        }
    }
    for(int i=0;i<3;i++)I.a[i][i]=1;
    gen[0]=1,gen[1]=root(p);
    gen[1]=Pow(gen[1],(p-1)/k);
    for(int i=2;i<k;i++)gen[i]=(1ll*gen[i-1]*gen[1])%p;
    for(int i=0;i<=2*k;i++)f[i]=omg(1ll*-i*(i-1)/2);
    for(int i=0;i<k;i++)g[i]=1ll*omg(1ll*i*(i-1)/2)*(Pow(W*omg(i)+I,l).a[x-1][y-1])%p;
    reverse(g,g+k+1);
    int lim=pre(2*k,k);
    Mul(lim,f,g);
    int inv=Pow(k,p-2);
    for(int t=0;t<k;t++){
        printf("%lld\n",1ll*omg(1ll*t*(t-1)/2)*inv%p*f[t+k]%p);
    }
    return 0;
}
/*
3 4 100 3 3 998244353
1 1 1
1 1 1
1 1 1

*/