题意
有一个((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} =
egin{cases}
1 i mod n = 0 \ 0 i mod n
eq 0 \
end{cases}
\这个等比数列求和即可证明；
\用这个来换掉式子中的取模：\
egin{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} imes 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_iomega_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)} imes 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 ",now); } return 0; }