HNOI2019 白兔之舞 dance

HNOI2019 白兔之舞 dance

---

显然(n=3)就是(n=1)的扩展版本，先来看看(n=1)怎么做。

令(W=w[1][1])，显然答案是：(ans_t=sum_{imod k=t}^{L}W^iinom{L}{i})

(=sum_{i=0}^{L}[k|(i-t)]W^iinom{L}{i})

这时用一个单位根反演。

回顾一下，单位根是fft时用到的东西，(omega_{n}=cosfrac{2pi}{n}+sinfrac{2pi}{n}i)，在膜(p)意义下，求出(p)的原根(g)，(omega_{n}=g^{frac{p-1}{n}})。

有一些单位根的性质：

(omega_n^i=omega_n^{imod n})

(omega_n^i=-omega_n^{i+frac{n}{2}})

单位根反演是这个东西：

(frac{1}{n}sum_{i=0}^{n-1}(omega_n^k)^i=[n|k])。

证明分类讨论：如果(n|k)，根据上面性质(omega_n^k=omega_n^0=1)；否则这是一个等比数列，公比为(omega_{n}^k)，求和为(frac{1-omega_n^{kn}}{1-omega_n^k})，上面的这个东西是(0)。

将单位根反演套进上面式子，代替([k|(i-t)])。

(=sum_{i=0}^{L}frac{1}{k}sum_{j=0}^{k-1}(omega_k^{i-t})^jW^iinom{L}{i})

(=frac{1}{k}sum_{j=0}^{k-1}omega_k^{-tj}sum_{i=0}^{L}W^iomega_k^{ij}inom{L}{i})

后面这个式子和t没什么关系了，可以统一计算，好像这是个二项式，可以直接算：

(=frac{1}{k}sum_{j=0}^{k-1}omega_k^{-tj}sum_{i=0}^{L}inom{L}{i}(Womega_k^j)^i1^{n-i})

(=frac{1}{k}sum_{j=0}^{k-1}omega_k^{-tj}(omega_k^jW+1)^L)

后面这东西只和j有关系，可以提前随便算出来，记为(c_j)

(=frac{1}{k}sum_{j=0}^{k-1}omega_k^{-tj}c_j)

然后神仙一波操作，(-tj=inom{t}{2}+inom{j}{2}-inom{t+j}{2})

(=frac{1}{k}sum_{j=0}^{k-1}omega_k^{inom{t}{2}+inom{j}{2}-inom{t+j}{2}}c_j)

(=frac{1}{k}omega_k^{inom{t}{2}}sum_{j=0}^{k-1}omega_k^{inom{j}{2}-inom{t+j}{2}}c_j)

(=frac{1}{k}omega_k^{inom{t}{2}}sum_{j=0}^{k-1}omega_k^{inom{j}{2}}c_jcdot omega_k^{-inom{t+j}{2}})

这就是个多项式乘法了，注意t+j那个多项式直接reverse一下，我真的越学越蠢了，不知道这个怎么搞自闭1h

然后因为模数不是ntt模数还要用MTT

好的(n=1)做完了，(n=3)实际上就是把(n=1)的(W)换成了输入的矩阵。于是只要修改一下(c_i)的求法，先算出((Begincdot(omega_k^jW)+I)^L)，Begin就是初始矩阵，只有(1,x)处是1，再取最终结果需要取的(1,y)处的值就行了

#include<bits/stdc++.h>
#define il inline
#define vd void
typedef long long ll;
il ll gi(){
    ll x=0,f=1;
    char ch=getchar();
    while(!isdigit(ch)){
        if(ch=='-')f=-1;
        ch=getchar();
    }
    while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
    return x*f;
}
const double pi=acos(-1);
int m,k,n,x,y,mod,G;
struct matrix{
    int s[3][3];matrix(){memset(s,0,sizeof s);}
};
matrix I;
int omega[65539];
il matrix operator+(const matrix&a,const matrix&b){
    matrix ret;
    for(int i=0;i<3;++i)
        for(int j=0;j<3;++j){
            ret.s[i][j]=a.s[i][j]+b.s[i][j];
            if(ret.s[i][j]>=mod)ret.s[i][j]-=mod;
        }
    return ret;
}
il matrix operator*(const matrix&a,const matrix&b){
    matrix ret;
    for(int j=0;j<3;++j)
        for(int i=0;i<3;++i)
            for(int k=0;k<3;++k)
                ret.s[i][k]=(ret.s[i][k]+1ll*a.s[i][j]*b.s[j][k])%mod;
    return ret;
}
il matrix operator*(const matrix&a,const int&b){
    matrix ret;
    for(int i=0;i<3;++i)
        for(int j=0;j<3;++j)
            ret.s[i][j]=1ll*a.s[i][j]*b%mod;
    return ret;
}
matrix w;
int A[262147],B[262147];
il int pow(int x,int y){
    int ret=1;
    while(y){
        if(y&1)ret=1ll*ret*x%mod;
        x=1ll*x*x%mod;y>>=1;
    }
    return ret;
}
il matrix pow(matrix x,int y){
    matrix ret=I;
    while(y){
        if(y&1)ret=ret*x;
        x=x*x;y>>=1;
    }
    return ret;
}
il int getrt(int x){
    static int p[50],o=0;
    for(int i=2,y=x-1;i<=y;++i)
        if(y%i==0){
            p[++o]=i;
            while(y%i==0)y/=i;
        }
    for(int g=2;;++g){
        bool yes=1;
        for(int j=1;j<=o;++j)if(pow(g,(mod-1)/p[j])==1){yes=0;break;}
        if(yes)return g;
    }
}
typedef std::complex<double> cp;
int rev[262147];cp omg[262147];
cp A1[262147],A2[262147],B1[262147],B2[262147];
il vd fft(cp*A,int n){
    for(int i=0;i<n;++i)if(rev[i]>i)std::swap(A[i],A[rev[i]]);
    for(int o=1;o<n;o<<=1)
        for(cp*p=A;p!=A+n;p+=o<<1)
            for(int i=0;i<o;++i){
                cp t=omg[n/(o<<1)*i]*p[i+o];
                p[i+o]=p[i]-t,p[i]+=t;
            }
}
int ans[262147];
int main(){
#ifdef XZZSB
    freopen("in.in","r",stdin);
    freopen("out.out","w",stdout);
#endif
    I.s[0][0]=I.s[1][1]=I.s[2][2]=1;
    m=gi(),k=gi(),n=gi(),x=gi()-1,y=gi()-1,mod=gi();
    omega[0]=1;omega[1]=pow(G=getrt(mod),(mod-1)/k);
    for(int i=2;i<k;++i)omega[i]=1ll*omega[1]*omega[i-1]%mod;
    for(int i=0;i<m;++i)for(int j=0;j<m;++j)w.s[i][j]=gi();
    for(int i=0;i<(k<<1|1);++i)A[i]=omega[(k-1ll*i*(i-1)/2%k)%k];
    matrix begin;begin.s[0][x]=1;
    for(int i=0;i<k;++i)B[i]=1ll*omega[1ll*i*(i-1)/2%k]*(begin*pow(w*omega[i]+I,n)).s[0][y]%mod;
    std::reverse(B,B+k+1);
    int N=1,lg=0;while(N<(k*3+5))N<<=1,++lg;
    for(int i=0;i<N;++i)omg[i]=(cp){cos(i*pi*2/N),sin(i*pi*2/N)};
    for(int i=0;i<N;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<lg-1);
    for(int i=0;i<N;++i)A1[i]=A[i]&32767,A2[i]=A[i]>>15;
    for(int i=0;i<N;++i)B1[i]=B[i]&32767,B2[i]=B[i]>>15;
    fft(A1,N),fft(B1,N);fft(A2,N),fft(B2,N);
    for(int i=0;i<N;++i){
        cp _A1=A1[i],_A2=A2[i],_B1=B1[i],_B2=B2[i];
        A1[i]=_A1*_B1;A2[i]=_A2*_B2;
        B1[i]=_A1*_B2;B2[i]=_A2*_B1;
    }
    for(int i=0;i<N;++i)omg[i]=(cp){cos(i*pi*2/N),-sin(i*pi*2/N)};
    fft(A1,N),fft(B1,N);fft(A2,N),fft(B2,N);
    for(int i=0;i<N;++i)A1[i]/=N;
    for(int i=0;i<N;++i)A2[i]/=N;
    for(int i=0;i<N;++i)B1[i]/=N;
    for(int i=0;i<N;++i)B2[i]/=N;
    for(int i=0;i<N;++i)A[i]=((ll)(A1[i].real()+0.5)%mod+1073741824ll*(((ll)(A2[i].real()+0.5))%mod)%mod+32768ll*(((ll)(B1[i].real()+0.5))%mod)%mod+32768ll*(((ll)(B2[i].real()+0.5))%mod)%mod)%mod;
    int invk=pow(k,mod-2);
    for(int i=0;i<k;++i)printf("%lld
",1ll*A[i+k]*invk%mod*omega[1ll*i*(i-1)/2%k]%mod);
    return 0;
}