题意
有一个\((L+1)*n\) 的网格图，初始时白兔在\((0,X)\) , 每次可以向横坐标递增，纵坐标随意的位置移动，两个位置之间的路径条数只取决于纵坐标，用\(w(i,j)\) 表示，如果要求白兔停下的点纵坐标为\(Y\) 依次输出移动的步数对\(k\) 取模为 $0 - k -1 $的方案数；
\(p\)为质数且$10^8 \lt p \lt 2^{30} , 1 \le n \le 3 , 1 \le x , y \le n , 0 \le w(i,j) \lt p , 1 \le k \le 65536 , \le L \le 10^8 $，
满足\(k\)是\(p-1\)的一个约数；
题解
即求：$ f_r = \sum_{i=0}{L}(L_i) (A^i)_{X,Y} [ i \ mod \ k = r] $
$Part \ 1 $
单位根反演：
\[根据求和引理：\\
\frac{1}{n} \sum_{j=0}^{n-1}\omega_n^{ij} =
\begin{cases}
1 \ \ \ i \ mod \ n = 0 \\ 0 \ \ \ i \ mod \ n \neq 0 \\
\end{cases}
\\这个等比数列求和即可证明；
\\用这个来换掉式子中的取模：\\
\begin{align}
&=\sum_{i=0}^{L}(^L_i)A^i_{x,y}\frac{1}{k}\sum_{j=0}^{k-1} \omega_k^{(i-r)j}\\
&=\frac{1}{k}\sum_{i=0}^{k-1}\omega_k^{-ir}\sum_{j=0}^{L}(^L_j)(A^j)_{x,y}\times \omega^{ij}_k\\
&=\frac{1}{k}\sum_{i=0}^{k-1}\omega_k^{-ir}(w^i_kA+I)^L_{x,y}\\
&=\frac{1}{k}\sum_{i=0}^{k-1}a_i\omega_k^{-ir}
\end{align}\]
\frac{1}{n} \sum_{j=0}^{n-1}\omega_n^{ij} =
\begin{cases}
1 \ \ \ i \ mod \ n = 0 \\ 0 \ \ \ i \ mod \ n \neq 0 \\
\end{cases}
\\这个等比数列求和即可证明；
\\用这个来换掉式子中的取模：\\
\begin{align}
&=\sum_{i=0}^{L}(^L_i)A^i_{x,y}\frac{1}{k}\sum_{j=0}^{k-1} \omega_k^{(i-r)j}\\
&=\frac{1}{k}\sum_{i=0}^{k-1}\omega_k^{-ir}\sum_{j=0}^{L}(^L_j)(A^j)_{x,y}\times \omega^{ij}_k\\
&=\frac{1}{k}\sum_{i=0}^{k-1}\omega_k^{-ir}(w^i_kA+I)^L_{x,y}\\
&=\frac{1}{k}\sum_{i=0}^{k-1}a_i\omega_k^{-ir}
\end{align}\]
当\(k = 2^x\)次方时，可以用裸的\(fft\)求出答案;
\(Part 2\)
但是没有必要，当 \(k\) 为任意数，由于\(k \ | \ mod - 1\)，所以\(\omega_k\)可以用原根的\((mod-1)/k\)来代替;
同时注意到： $ ab = (^{a+b}_2) - (^a_2) - (^b_2) $
$ f_r = \frac{1}{k}w_k{(r_2)} \sum_{i=0}^{k-1}a_i \omega _k{(i_2)} \times \omega_k{-({i+r}_2)} $
-
直接reverse+mtt即可;
//注意mtt的精度;
#include<bits/stdc++.h>
#define ll long long
#define ld double
using namespace std;
const int N=1<<20;
int n,k,l,X,Y,mod,a[N],b[N],A[N],B[N],C[N],fac[N],tot;
char gc(){
static char*p1,*p2,s[1000000];
if(p1==p2)p2=(p1=s)+fread(s,1,1000000,stdin);
return(p1==p2)?EOF:*p1++;
}
int rd(){
int x=0;char c=gc();
while(c<'0'||c>'9')c=gc();
while(c>='0'&&c<='9')x=(x<<1)+(x<<3)+c-'0',c=gc();
return x;
}
int pw(int x,int y){
int re=1;
while(y){
if(y&1)re=(ll)re*x%mod;
y>>=1;x=(ll)x*x%mod;
}
return re;
}
int find_rt(int p){
int tmp=p-1;
for(int i=2;i*i<=tmp;++i)if(tmp%i==0){
fac[++tot]=i;
while(tmp%i==0)tmp/=i;
}
if(tmp!=1)fac[++tot]=tmp;
for(int i=2;;++i){
int fg=0;
for(int j=1;j<=tot;++j)if(pw(i,(p-1)/fac[j])==1){fg=1;break;}
if(!fg)return i;
}
return 0;
}
struct Mat{
int a[3][3];
void init(){
memset(a,0,sizeof(a));
for(int i=0;i<n;++i)a[i][i]=1;
}
Mat operator *(const int&A){
Mat re;
for(int i=0;i<n;++i)
for(int j=0;j<n;++j)re.a[i][j]=(ll)a[i][j]*A%mod;
return re;
}
Mat operator *(const Mat&A){
Mat re;
for(int i=0;i<n;++i)
for(int j=0;j<n;++j){
re.a[i][j]=0;
for(int k=0;k<n;++k){
re.a[i][j]+=(ll)a[i][k]*A.a[k][j]%mod;
if(re.a[i][j]>=mod)re.a[i][j]-=mod;
}
}
return re;
}
Mat operator +(const Mat&A){
Mat re;
for(int i=0;i<n;++i)
for(int j=0;j<n;++j){
re.a[i][j]=a[i][j]+A.a[i][j];
if(re.a[i][j]>=mod)re.a[i][j]-=mod;
}
return re;
}
}mp,I;
int pw(Mat x,int y){
Mat re;re.init();
while(y){
if(y&1)re=re*x;
y>>=1;x=x*x;
}
return re.a[X][Y];
}
struct com{
ld x,y;
com(ld _x=0,ld _y=0):x(_x),y(_y){};
com operator +(const com&A)const{return com(x+A.x,y+A.y);}
com operator -(const com&A)const{return com(x-A.x,y-A.y);}
com operator *(const com&A)const{return com(x*A.x-y*A.y,x*A.y+y*A.x);}
com operator /(const ld&A)const{return com(x/A,y/A);}
com operator !()const{return com(x,-y);}
};
namespace poly{
const ld pi=acos(-1);
int L,len,rev[N];com Wn[N];
void init(int l){
for(L=0,len=1;len<=l<<1;len<<=1,++L);
for(int i=0;i<len;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
//Wn[0]=com(1,0);Wn[1]=com(cos(pi/len),sin(pi/len));
//for(int i=2;i<len;++i)Wn[i]=Wn[i-1]*Wn[1];
for(int i=0;i<len;++i)Wn[i]=com(cos(pi*i/len),sin(pi*i/len));
}
void fft(com*A,int f){
for(int i=0;i<len;++i)if(i<rev[i])swap(A[i],A[rev[i]]);
for(int i=1;i<len;i<<=1){
for(int j=0;j<len;j+=i<<1){
for(int k=0;k<i;++k){
com w=~f?Wn[len/i*k]:!Wn[len/i*k];//
com x=A[j+k],y=w*A[j+k+i];
A[j+k]=x+y;A[j+k+i]=x-y;
}
}
}
if(!~f)for(int i=0;i<len;++i)A[i]=A[i]/len;
}
void mtt(int*A,int*B,int*C){
static com t1[N],t2[N],t3[N],t4[N];
int all=(1<<15)-1;
for(int i=0;i<len;++i){
t1[i]=com(A[i]>>15,A[i]&all);
t2[i]=com(B[i]>>15,B[i]&all);
}
fft(t1,1);fft(t2,1);
for(int i=0;i<len;++i){
com x1=t1[i],y1=!t1[len-i&len-1];
com x2=t2[i],y2=!t2[len-i&len-1];
t3[i] = (x1+y1)*x2*com(0.5,0);//(x1+y1)/2(x2+y2)/2 + (x1+y1)/2(x2-y2)/2i * i ;
t4[i] = (x1-y1)*x2*com(0,-0.5);//(x1-y1)/2i(x2+y2)/2 + (x1-y1)/2i(x2-y2)/2i * i ;
}
fft(t3,-1);fft(t4,-1);
for(int i=0;i<len;++i){
C[i] = (((ll)(t3[i].x+0.5)%mod <<30)%mod +
((ll)(t3[i].y+0.5)%mod <<15)%mod +
((ll)(t4[i].x+0.5)%mod <<15)%mod +
(ll)(t4[i].y+0.5)%mod ) %mod ;
}
}
}
int main(){
// freopen("dance.in","r",stdin);
// freopen("dance.out","w",stdout); n=rd();k=rd();l=rd();X=rd()-1;Y=rd()-1;mod=rd();I.init();
for(int i=0;i<n;++i)for(int j=0;j<n;++j)mp.a[i][j]=rd();
int G=find_rt(mod),wk=pw(G,(mod-1)/k); for(int i=0,w=1;i<k;++i,w=(ll)w*wk%mod)a[i]=pw((mp*w+I),l);
for(int i=0;i<=k-1<<1;++i)b[i]=pw(wk,(ll)i*(i-1)/2%(mod-1));
for(int i=0;i<k;++i)A[i]=(ll)a[i]*b[i]%mod;
for(int i=0;i<=k-1<<1;++i)B[i]=pw(b[(k-1<<1)-i],mod-2); poly::init(k-1<<1);
poly::mtt(A,B,C);
int iv=pw(k,mod-2);
for(int i=0;i<k;++i){
int now = (ll)iv*b[i]%mod*C[(k-1<<1)-i]%mod;
printf("%d\n",now);
} return 0;
}