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;
}