HNOI2019 白兔之舞 dance
显然\(n=3\)就是\(n=1\)的扩展版本，先来看看\(n=1\)怎么做。
令\(W=w[1][1]\)，显然答案是：\(ans_t=\sum_{i\mod k=t}^{L}W^i\binom{L}{i}\)
\(=\sum_{i=0}^{L}[k|(i-t)]W^i\binom{L}{i}\)
这时用一个单位根反演。
回顾一下，单位根是fft时用到的东西，\(\omega_{n}=\cos\frac{2\pi}{n}+\sin\frac{2\pi}{n}i\)，在膜\(p\)意义下，求出\(p\)的原根\(g\)，\(\omega_{n}=g^{\frac{p-1}{n}}\)。
有一些单位根的性质：
\(\omega_n^i=\omega_n^{i\mod 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^i\binom{L}{i}\)
\(=\frac{1}{k}\sum_{j=0}^{k-1}\omega_k^{-tj}\sum_{i=0}^{L}W^i\omega_k^{ij}\binom{L}{i}\)
后面这个式子和t没什么关系了，可以统一计算，好像这是个二项式，可以直接算：
\(=\frac{1}{k}\sum_{j=0}^{k-1}\omega_k^{-tj}\sum_{i=0}^{L}\binom{L}{i}(W\omega_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=\binom{t}{2}+\binom{j}{2}-\binom{t+j}{2}\)
\(=\frac{1}{k}\sum_{j=0}^{k-1}\omega_k^{\binom{t}{2}+\binom{j}{2}-\binom{t+j}{2}}c_j\)
\(=\frac{1}{k}\omega_k^{\binom{t}{2}}\sum_{j=0}^{k-1}\omega_k^{\binom{j}{2}-\binom{t+j}{2}}c_j\)
\(=\frac{1}{k}\omega_k^{\binom{t}{2}}\sum_{j=0}^{k-1}\omega_k^{\binom{j}{2}}c_j\cdot \omega_k^{-\binom{t+j}{2}}\)
这就是个多项式乘法了，注意t+j那个多项式直接reverse一下，我真的越学越蠢了，不知道这个怎么搞自闭1h
然后因为模数不是ntt模数还要用MTT
好的\(n=1\)做完了，\(n=3\)实际上就是把\(n=1\)的\(W\)换成了输入的矩阵。于是只要修改一下\(c_i\)的求法，先算出\((Begin\cdot(\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\n",1ll*A[i+k]*invk%mod*omega[1ll*i*(i-1)/2%k]%mod);
return 0;
}