【HNOI2019】白兔之舞【组合数学】【矩阵快速幂】【单位根反演】【Chirp Z-Transform】【原根】【MTT】

题意：有一张 $(L+1)×n(L+1)\times n$ 个点的有向图，每个结点有二元组 $(x,y)(0≤x≤L,1≤y≤n)(x,y)~(0\leq x\leq L,1\leq y\leq n)$ 表示。对于所有 $u_1,v_1),(u_2,v_2)$ ，若 $u_1<u_2$ ,则前者向后者连 $w_{v_1,v_2}$ 条边。对所有 $0≤t<k0\leq t<k$ ，你需要求出从 $(0, x)$ 走模 $k$ 余 $t$ 步到达任意一个第二维为 $y$ 的点的方案数。答案对 $p$ 取模。

$L≤108,n≤3,k≤216,108<p<230,p∈prime,k∣p−1L\leq 10^8,n\leq 3,k\leq 2^{16},10^8<p<2^{30},p\in prime,k\mid p-1$

神奇的题

显然是单位根反演，所以先想办法搞出走到位置 $i$ 的生成函数。

显然是个矩阵的形式。看成只有 $n$ 个点在上面跑，邻接矩阵为 $M$

那么走一步的生成函数为

$s(x)=∑i=1∞Mxis(x)=\sum_{i=1}^{\infin}Mx^i$

答案的生成函数为

$f(x)=∑i=0∞(∑j=0L[xj]si(x))xif(x)=\sum_{i=0}^{\infin}\left(\sum_{j=0}^L[x^j]s^i(x)\right)x^i$

看上去很不可做，但实际上中间那坨看起来形式很简单。冷静分析，实际上 $M$ 的次数就是 $i$ ，前面的系数就是选 $i$ 个正整数和不超过 $L$ ,等价于选 $i + 1$ 个正整数不超过 $L + 1$ ，即

$f(x)=∑i=0L(Li)Mixi=(Mx+1)Lf(x)=\sum_{i=0}^L\binom LiM^ix^i\\=(Mx+1)^L$

然后就可以单位根反演了

$anst=1k∑i=0k−1ωk−it(Mωki+1)Lans_t=\frac 1k\sum_{i=0}^{k-1}\omega_k^{-it}(M\omega_{k}^i+1)^L$

我们发现我们只求这个东西的 $x$ 行 $y$ 列，所以可以设 $ai=[(Mωki+1)L]x,ya_i=\left[(M\omega_k^i+1)^L\right]_{x,y}$

然后

$anst=1k∑i=0k−1ωk−itaians_t=\frac 1k\sum_{i=0}^{k-1}\omega_k^{-it}a_i$

这是个 CZT 的形式，拆开就可以了

$anst=1k∑i=0k−1ωk(i2)+(t2)−(i+t2)aians_t=\frac 1k\sum_{i=0}^{k-1}\omega_k^{\binom i2+\binom t2-\binom{i+t}2}a_i$

$anst=ωk(t2)k∑i=0k−1ωk(i2)aiωk−(i+t2)ans_t=\dfrac{\omega_k^{\binom t2}}{k}\sum_{i=0}^{k-1}\omega_k^{\binom i2}a_i\omega_k^{-\binom{i+t}2}$

直接做减法卷积就可以了。

求个原根出来求单位元，然后 MTT 即可。

因为要写的东西太多了，很多地方都写的很浪，仅供参考。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cctype>
#include <cmath>
#include <algorithm>
#define MAXN (1<<18)+5
#define double long double
using namespace std;
typedef long long ll;
inline int read()
{int ans=0;char c=getchar();while (!isdigit(c)) c=getchar();while (isdigit(c)) ans=(ans<<3)+(ans<<1)+(c^48),c=getchar();return ans;
}
int n,k,L,x,y,M,P,g;
inline int qpow(int a,int p)
{int ans=1;while (p){if (p&1) ans=(ll)ans*a%P;a=(ll)a*a%P,p>>=1;}return ans;
}
const double Pi=acos(-1.0);
struct complex{double x,y;inline complex(const double& x=0,const double& y=0):x(x),y(y){}};
inline complex operator +(const complex& a,const complex& b){return complex(a.x+b.x,a.y+b.y);}
inline complex operator -(const complex& a,const complex& b){return complex(a.x-b.x,a.y-b.y);}
inline complex operator *(const complex& a,const complex& b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
inline complex adj(const complex& a){return complex(a.x,-a.y);}
int l,lim,r[MAXN];
inline void init(){lim=1<<l;for (int i=0;i<lim;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));}
complex rt[2][24];
void fft(complex* a,int type)
{for (int i=0;i<lim;i++) if (i<r[i]) swap(a[i],a[r[i]]);for (int L=0;L<l;L++){int mid=1<<L,len=mid<<1;complex Wn=rt[type][L+1];for (int s=0;s<lim;s+=len){complex w(1,0);for (int k=0;k<mid;k++,w=w*Wn){complex x=a[s+k],y=w*a[s+mid+k];a[s+k]=x+y,a[s+mid+k]=x-y;}}}if (type) for (int i=0;i<lim;i++) a[i].x/=lim,a[i].y/=lim;
}
complex A[MAXN],B[MAXN],C[MAXN],D[MAXN];
void fmul(int* F,int* G,int n,int m)
{M=sqrt(P);for (l=0;(1<<l)<n+m;l++);init();for (int i=0;i<n;i++) A[i]=complex(F[i]/M,F[i]%M);for (int i=0;i<m;i++) C[i]=complex(G[i]/M,G[i]%M);fft(A,0),fft(C,0);B[0]=adj(A[0]),D[0]=adj(C[0]);for (int i=1;i<lim;i++) B[i]=adj(A[lim-i]),D[i]=adj(C[lim-i]);for (int i=0;i<lim;i++){complex a=A[i],b=B[i],c=C[i],d=D[i];A[i]=(a+b)*complex(0.5,0);B[i]=(a-b)*complex(0,-0.5);C[i]=(c+d)*complex(0.5,0);D[i]=(c-d)*complex(0,-0.5);}for (int i=0;i<lim;i++){complex a=A[i],b=B[i],c=C[i],d=D[i];A[i]=a*c+a*d*complex(0,1);B[i]=b*c+b*d*complex(0,1);}fft(A,1),fft(B,1);for (int i=0;i<n+m-1;i++){ll a=A[i].x+0.5,b=A[i].y+B[i].x+0.5,c=B[i].y+0.5;a%=P,b%=P,c%=P;F[i]=(a*M%P*M%P+b*M%P+c)%P;}
}
struct mat
{int e[3][3];inline mat(const int& v=0){memset(e,0,sizeof(e));for (int i=0;i<3;i++) e[i][i]=v;	}	inline int* operator [](const int& i){return e[i];}inline const int* operator [](const int& i)const{return e[i];}
}W;
inline mat operator +(const mat& a,const mat& b)
{mat c;for (int i=0;i<3;i++)for (int j=0;j<3;j++)c[i][j]=(a[i][j]+b[i][j])%P;return c;
}
inline mat operator *(const mat& a,const mat& b)
{mat c;for (int i=0;i<3;i++)for (int k=0;k<3;k++)for (int j=0;j<3;j++)c[i][j]=(c[i][j]+(ll)a[i][k]*b[k][j])%P;return c;
}
inline mat qpow(mat a,int p)
{mat ans(1);while (p){if (p&1) ans=ans*a;a=a*a,p>>=1;}return ans;
}
inline bool check(int g)
{int x=P-1;for (int i=2;i*i<=x;i++)if (x%i==0){if (qpow(g,(P-1)/i)==1) return false;while (x%i==0) x/=i;}if (x>1&&qpow(g,(P-1)/x)==1) return false; return true;
}
inline int sum(int x){return (ll)x*(x-1)/2%k;}
int F[MAXN],G[MAXN];
int main()
{for (int i=0;i<20;i++){double a=2*Pi/(1<<i);rt[1][i]=adj(rt[0][i]=complex(cos(a),sin(a)));}n=read(),k=read(),L=read(),x=read()-1,y=read()-1,P=read();while (!check(++g));for (int i=0;i<n;i++) for (int j=0;j<n;j++) W[i][j]=read();int Wn=qpow(g,(P-1)/k);for (int i=0,w=1;i<k;i++,w=(ll)w*Wn%P) F[i]=(ll)qpow(W*mat(w)+mat(1),L)[x][y]*qpow(Wn,sum(i))%P;for (int i=0;i<2*k;i++)G[i]=qpow(Wn,k-sum(i));
//	for (int t=0;t<k;t++)
//	{
//		int ans=0;
//		for (int i=0;i<k;i++) ans=(ans+(ll)F[i]*G[i+t])%P;
//		printf("%d ",ans);
//	}
//	puts("");reverse(F,F+k);	fmul(F,G,k,2*k);
//	for (int i=k-1;i<2*k-1;i++) printf("%d ",F[i]);
//	puts("");int t=qpow(k,P-2);for (int i=k-1;i<2*k-1;i++) printf("%lld\n",(ll)F[i]*qpow(Wn,sum(i-k+1))%P*t%P);return 0;
}