单位根反演不会啊怎么搞 $FFT$ 吧，还是了解了单位根反演后才可以搞的好吧…… 居然有人吐槽我说我学了 $FFT$ 但是不会运用？！，嘤嘤嘤打击有些大……
实际上所谓的单位根反演就是这个东西：
$$\frac{1}{n}\sum_{i=0}^{n-1}(\omega_n^d)^i=[n|d]$$
回到题目，我们先考虑正解的简化版—— $n=1$ 的版本，我们先定义 $W=w[1][1]$ 。
现在对于每一个 $t$ 的答案显然为 $\sum_{i=0}^{L}[i\% k=t] W^i (^L_i)$
这个式子显然等于 $\sum_{i=0}^{L}[k|(i-t)] w^i (^L_i)$ 。会发现 $[k|(i-t)]$ 和上面单位根反演的 $[n|d]$ 一样，于是我们尝试将单位根反演的式子带进去。
$$
=\sum_{i=0}^{L}\frac{1}{k}\sum_{j=0}^{k-1}(\omega_k^{i-t})^j W^i \binom{L}{i}\\
=\frac{1}{k}\sum_{j=0}^{k-1}\omega_k^{-tj}\sum_{i=0}^{L} \omega_k^{ij} W^i \binom{L}{i}\\
=\frac{1}{k}\sum_{j=0}^{k-1}\omega_k^{-tj}\sum_{i=0}^{L} \binom{L}{i}(\omega_k^{j} W)^i\\
=\frac{1}{k}\sum_{j=0}^{k-1}\omega_k^{-tj}\sum_{i=0}^{L} \binom{L}{i}(\omega_k^{j} W)^i 1^{n-i}\\
=\frac{1}{k}\sum_{j=0}^{k-1}\omega_k^{-tj}(\omega_k^{j} W+1)^L
$$
后面的 $(\omega_k^{j} W+1)^L$ 显然可以预处理，记为 $num_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}}num_j\\
=\frac{1}{k}\omega_k^{\binom{t}{2}}\sum_{j=0}^{k-1}num_j\omega_k^{\binom{j}{2}}\cdot\omega_k^{-\binom{t+j}{2}}$$
后面的式子可以用 $FFT$ 加速，但是值域太大这里需要用到 $MTT$ 。现在就有 $40$ 分了，接下来考虑 $n>1$ 的情况。
我们建矩阵，然后会发现 $n>1$ 仅会对 $num_j$ 的计算方式有变化。
我们定义一个 $begin$ 矩阵，该矩阵只有 $(0,x)$ 位置上有值且值为 $1$ ，也就是说这是白兔的起点。那么最后我们需要留下来的也就是矩阵的 $(0,y)$ ，因为只有在第二维为 $y$ 是才会计入答案。
嗯，差不多可以这样写：
begin.c[0][x]=1;
for(int i=0;i<k;++i)
num[i]=(begin*pow(w*num[i]+I,n)).c[0][y]%MOD;
/*w 就是上文中的 W，不过这里是矩阵*/
/*I 是矩阵中的单位'1'*/
Code:
#include <cmath>
#include <cstdio>
#include <string>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
typedef long long ll;
const int N=65536;
const double PI=acos(-1);
int m,k,n,x,y,MOD,G,num[N],A[N<<2],B[N<<2];
namespace OI {
#define F(x,i,j) for((x)=(i);(x)<=(j);++(x))
inline ll IN() {
char ch;bool flag=0;ll x=0;
while(ch=getchar(),!isdigit(ch)) if(ch=='-') flag=1;
while(isdigit(ch)) x=x*10+ch-'0',ch=getchar();
if(flag) x=-x;return x;
}
struct matrix {int c[3][3];matrix(){memset(c,0,sizeof(c));}};
matrix operator + (const matrix&a,const matrix&b) {
matrix ans;int i,j;F(i,0,2)F(j,0,2) {
ans.c[i][j]=a.c[i][j]+b.c[i][j];
if(ans.c[i][j]>=MOD) ans.c[i][j]-=MOD;
}return ans;
}
matrix operator * (const matrix&a,const matrix&b) {
matrix ans;int i,j,k;F(i,0,2)F(j,0,2)F(k,0,2)
ans.c[i][k]=(ans.c[i][k]+1ll*a.c[i][j]*b.c[j][k])%MOD;
return ans;
}
matrix operator * (const matrix&a,const int&b) {
matrix ans;int i,j;F(i,0,2)F(j,0,2)ans.c[i][j]=1ll*a.c[i][j]*b%MOD;
return ans;
}
struct complex{complex(long double a=0,long double b=0){x=a,y=b;}long double x,y;};
complex operator + (complex a,complex b){return complex(a.x+b.x,a.y+b.y);}
complex operator - (complex a,complex b){return complex(a.x-b.x,a.y-b.y);}
complex operator * (complex a,complex b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
matrix I;
inline int pow(int x,int y) {
int res=1;for(;y;y>>=1,x=1ll*x*x%MOD) if(y&1) res=1ll*res*x%MOD;
return res%MOD;
}
inline matrix pow(matrix x,int y) {
matrix res=I;for(;y;y>>=1,x=x*x) if(y&1) res=res*x;
return res;
}
}using namespace OI;
namespace MTT {
#define BLOCK 32768
int limit=1,cnt=0,filp[N<<2],Ans[N<<2];
complex A1[N<<2],B1[N<<2],A2[N<<2],B2[N<<2],X[N<<2],omg[N<<2];
inline void fft(complex *f,short inv){
for(int i=0;i<limit;++i)if(i<filp[i])std::swap(f[i],f[filp[i]]);
for(int p=1;p<limit;p<<=1)
for(complex *a=f;a!=f+limit;a+=(p<<1))
for(int l=0;l<p;++l){
complex t=a[l+p]*omg[limit/(p<<1)*l];
a[l+p]=a[l]-t,a[l]=a[l]+t;
}
}
inline void mtt(int *A,int *B){
while(limit<(k*3+5)) limit<<=1,++cnt;
for(int i=0;i<limit;++i) filp[i]=(filp[i>>1]>>1)|((i&1)<<(cnt-1));
for(int i=0;i<limit;++i) A1[i].x=A[i]&(BLOCK-1),A2[i].x=A[i]>>15;
for(int i=0;i<limit;++i) B1[i].x=B[i]&(BLOCK-1),B2[i].x=B[i]>>15;
for(int i=0;i<limit;++i) omg[i]=(complex){cos(i*PI*2/limit),sin(i*PI*2/limit)};
fft(A1,1),fft(B1,1);fft(A2,1),fft(B2,1);
for(int i=0;i<limit;++i){
complex 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<limit;++i) omg[i]=(complex){cos(i*PI*2/limit),-sin(i*PI*2/limit)};
fft(A1,-1),fft(B1,-1);fft(A2,-1),fft(B2,-1);
for(int i=0;i<limit;++i)
A1[i].x/=limit,A2[i].x/=limit,B1[i].x/=limit,B2[i].x/=limit;
for(int i=0;i<limit;++i)
Ans[i]=((ll)(A1[i].x+0.5)%MOD+1073741824ll*((ll)(A2[i].x+0.5)%MOD)%MOD+
32768ll*((ll)(B1[i].x+0.5)%MOD)%MOD+32768ll*((ll)(B2[i].x+0.5)%MOD)%MOD)%MOD;
}
}using namespace MTT;
int divisor[105],tot;
inline int get_G() {/*获取原根*/
for(int i=2,S=MOD-1;i<=S;++i)
if(S%i==0) {divisor[++tot]=i;while(!(S%i)) S/=i;}
for(int g=2;;++g) {
bool ok=true;
for(int j=1;j<=tot;++j) if(pow(g,(MOD-1)/divisor[j])==1) {ok=false;break;}
if(ok) return g;
}
}
matrix w,s;
int main() {
I.c[0][0]=I.c[1][1]=I.c[2][2]=1;
m=IN(),k=IN(),n=IN(),x=IN(),y=IN(),MOD=IN();--x,--y;
/*num 其实就是上文中的单位根，这里预处理一下计算方便些*/
num[0]=1,num[1]=pow(G=get_G(),(MOD-1)/k);
for(int i=2;i<k;++i) num[i]=1ll*num[1]*num[i-1]%MOD;
for(int i=0;i<m;++i) for(int j=0;j<m;++j) w.c[i][j]=IN();
for(int i=0;i<(k<<1|1);++i) A[i]=num[(k-1ll*i*(i-1)/2%k)%k];
s.c[0][x]=1;
for(int i=0;i<k;++i) B[i]=1ll*num[1ll*i*(i-1)/2%k]*(s*pow(w*num[i]+I,n)).c[0][y]%MOD;
/*计算后面两个多项式的值*/
reverse(B,B+k+1),mtt(A,B);
int invk=pow(k,MOD-2);
for(int i=0;i<k;++i) printf("%lld\n",1ll*Ans[i+k]*invk%MOD*num[1ll*i*(i-1)/2%k]%MOD);
/*计算答案*/
return 0;
}
0 条评论