Loj 3058. 「HNOI2019」白兔之舞
题目描述
有一张顶点数为 ((L+1) imes n) 的有向图。这张图的每个顶点由一个二元组 ((u,v)) 表示 ((0le ule L,1le vle n))。这张图不是简单图，对于任意两个顶点 ((u_1,v_1),(u_2,v_2))，如果 (u_1<u_2)，则从 ((u_1,v_1)) 到 ((u_2,v_2)) 一共有 (w(v_1,v_2)) 条不同的边，如果 (u_1ge u_2) 则没有边。
白兔将在这张图上上演一支舞曲。白兔初始时位于该有向图的顶点 ((0,x))。
白兔将会跳若干步。每一步，白兔会从当前顶点沿任意一条出边跳到下一个顶点。白兔可以在任意时候停止跳舞（也可以没有跳就直接结束）。当到达第一维为 (L) 的顶点就不得不停止，因为该顶点没有出边。
假设白兔停止时，跳了 (m) 步，白兔会把这只舞曲给记录下来成为一个序列。序列的第 (i) 个元素为它第 (i) 步经过的边。
问题来了：给定正整数 (k) 和 (y (1le yle n))，对于每个 (t (0le t<k))，求有多少种舞曲（假设其长度为 (m)）满足 (m mod k=t)，且白兔最后停在了坐标第二维为 (y) 的顶点？
两支舞曲不同定义为它们的长度（(m)）不同或者存在某一步它们所走的边不同。
输出的结果对 (p) 取模。
输入格式
第一行六个用空格隔开的整数 (n,k,L,x,y,p)。
接下来 (n) 行，每行有 (n) 个用空格隔开的整数，第 (i) 行的第 (j) 个数表示 (w(i,j))。
输出格式
依次输出 (k) 行，每行一个数表示答案对 (p) 取模的结果。
样例
样例输入 1
2 2 3 1 1 998244353
2 1
1 0
样例输出 1
16
18
样例输入 2
3 4 100 1 3 998244353
1 1 1
1 1 1
1 1 1
样例输出 2
506551216
528858062
469849094
818871911
数据范围与提示
对于全部数据，(p) 为一个质数，(10^8<p<2^{30})，(1le nle 3)，(1le xle n)，(1le yle n)，(0le w(i,j)<p)，(1le kle 65536)，(k) 为 (p-1) 的约数，(1le Lle 10^8)。
对于每组测试点，特殊限制如下：
- 测试点 (1,2)：(Lle 10^5)；
- 测试点 (3)：(n=1,w(1,1)=1)，(k) 的最大质因子为 (2)；
- 测试点 (4)：(n=1)，(k) 的最大质因子为 (2)；
- 测试点 (5)：(n=1,w(1,1)=1)；
- 测试点 (6)：(n=1)；
- 测试点 (7,8)：(k) 的最大质因子为 (2)。
(\)
湖南的题好神啊，只能(Orz)了。
先考虑从(x)走(L)步到(y)，假设每次向右只走一格的方案数。很容易写出一个(DP)方程，然后用矩阵快速幂优化。我们读进来的矩阵(A)就是转移矩阵。所以方案数就是([A^i]_{x,y})。
计算答案的时候就枚举走了(i)步：
[ans_t=sum_{i=0}^L[i\%k==t]inom{L}{i}[A^i]_{x,y}\
=sum_{i=0}^Lfrac{1}{k}sum_{j=0}^{k-1}omega_k^{(i-t)*j}inom{L}{i}[A^i]_{x,y}\
=frac{1}{k}sum_{j=0}^{k-1}omega_k^{-jt}
sum_{i=0}^Linom{L}{i}(omega_k^{ij} [A^i]_{x,y})\
=frac{1}{k}sum_{j=0}^{k-1}omega_k^{-jt}[(omega_k^jA+I)^L]_{x,y}\
]
其中：
[I=egin{bmatrix}
1&0&0\
0&1&0\
0&0&1\
end{bmatrix}
]
上面最后一步变换实质上就是二项式定理(sum_{i=0}^ninom{n}{i}x_i=(1+x)^n)的运用。
设：
[f_j=[(omega_k^jA+I)^L]_{x,y}
]
则：
[ans_t=frac{1}{k}sum_{j=0}^{k-1}omega_{k}^{-jt}f_j\
]
然后就是一个很牛逼的套路：
[-tj=inom{t}{2}+inom{j}{2}-inom{t+j}{2}
]
所以：
[ans_t=frac{1}{k}sum_{j=0}^{k-1}omega_k^{inom{t}{2}+inom{j}{2}-inom{t+j}{2}}f_j\
=frac{1}{k}omega_k^{inom{t}{2}}sum_{j=0}^{k-1}omega_k^{inom{j}{2}}f_jcdotomega_k^{-inom{t+j}{2}}
]
(ans_i)的生成函数可以写成下面两个函数卷积。
[F(n)=omega_k^{inom{i}{2}}cdot f_n\
G(n)=omega_k^{-inom{2*k-1-n}{2}}\
]
代码：
#include<bits/stdc++.h>
#define ll long long
#define K 200005
using namespace std;
inline ll Get() {ll x=0,f=1;char ch=getchar();while(ch<'0'||ch>'9') {if(ch=='-') f=-1;ch=getchar();}while('0'<=ch&&ch<='9') {x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}return x*f;}
ll n,k,L,x,y,mod;
ll g,W,w[K];
ll f[K];
ll ksm(ll t,ll x) {
ll ans=1;
for(;x;x>>=1,t=t*t%mod)
if(x&1) ans=ans*t%mod;
return ans;
}
vector<int>fac;
void work(int t) {
for(int i=2,maxx=sqrt(t);i<=maxx;i++) {
if(t%i==0) fac.push_back(i);
while(t%i==0) t/=i;
}
if(t>1) fac.push_back(t);
}
ll Find_root() {
for(int i=2;i<mod;i++) {
if(ksm(i,mod-1)!=1) continue ;
int flag=1;
for(int j=0;j<fac.size();j++) {
if(ksm(i,(mod-1)/fac[j])==1) {
flag=0;
break;
}
}
if(flag) return i;
}
}
struct matrix {
ll a[4][4];
void Init() {memset(a,0,sizeof(a));}
}A,B;
matrix operator *(const matrix &a,const matrix &b) {
static matrix tem;
tem.Init();
for(int i=1;i<=n;i++) {
for(int j=1;j<=n;j++) {
for(int k=1;k<=n;k++)
tem.a[i][j]+=a.a[i][k]*b.a[k][j]%mod;
tem.a[i][j]%=mod;
}
}
return tem;
}
matrix ksm(matrix g,ll x) {
static matrix ans;
ans.Init();
for(int i=1;i<=n;i++) ans.a[i][i]=1;
for(;x;x>>=1,g=g*g)
if(x&1) ans=ans*g;
return ans;
}
ll C_2(ll n) {return n*(n-1)/2%k;}
struct Com {
long double r,v;
Com() {}
Com(long double _r,long double _v) {r=_r,v=_v;}
};
Com operator *(const Com &a,const Com &b) {return Com(a.r*b.r-a.v*b.v,a.r*b.v+a.v*b.r);}
Com operator +(const Com &a,const Com &b) {return Com(a.r+b.r,a.v+b.v);}
Com operator -(const Com &a,const Com &b) {return Com(a.r-b.r,a.v-b.v);}
Com operator /(const Com &a,const long double &b) {return Com(a.r/b,a.v/b);}
const long double pi=acos(-1);
void FFT(Com *a,int d,int flag) {
int n=1<<d;
static int rev[K<<2];
for(int i=0;i<n;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<d-1);
for(int i=0;i<n;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int s=1;s<=d;s++) {
int len=1<<s,mid=len>>1;
Com w(cos(2*flag*pi/len),sin(2*flag*pi/len));
for(int i=0;i<n;i+=len) {
Com t(1,0);
for(int j=0;j<mid;j++,t=t*w) {
Com u=a[i+j],v=a[i+j+mid]*t;
a[i+j]=u+v;
a[i+j+mid]=u-v;
}
}
}
if(flag==-1) {for(int i=0;i<n;i++) a[i]=a[i]/n;}
}
ll F[K<<2],G[K<<2],ans[K<<2];
void mul(ll *a,ll *b,ll *c,int d) {
Com A[K<<2],B[K<<2],C[K<<2],D[K<<2];
Com f1[K<<2],f2[K<<2],f3[K<<2],f4[K<<2];
ll maxx=sqrt(mod);
int n=1<<d;
for(int i=0;i<n;i++) {
A[i]=Com(a[i]/maxx,0);
B[i]=Com(a[i]%maxx,0);
C[i]=Com(b[i]/maxx,0);
D[i]=Com(b[i]%maxx,0);
}
FFT(A,d,1),FFT(B,d,1);
FFT(C,d,1),FFT(D,d,1);
for(int i=0;i<n;i++) {
f1[i]=A[i]*C[i];
f2[i]=A[i]*D[i];
f3[i]=B[i]*C[i];
f4[i]=B[i]*D[i];
}
FFT(f1,d,-1),FFT(f2,d,-1);
FFT(f3,d,-1),FFT(f4,d,-1);
for(int i=0;i<n;i++) {
c[i]=((ll)(f1[i].r+0.5)*maxx%mod*maxx%mod+(ll)(f2[i].r+0.5)*maxx%mod+(ll)(f3[i].r+0.5)*maxx+(ll)(f4[i].r+0.5))%mod;
}
}
int main() {
n=Get(),k=Get(),L=Get(),x=Get(),y=Get(),mod=Get();
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
A.a[i][j]=Get();
work(mod-1);
g=Find_root();
w[0]=1;
w[1]=ksm(g,(mod-1)/k);
for(int i=2;i<k;i++) w[i]=w[i-1]*w[1]%mod;
for(int j=0;j<k;j++) {
B=A;
for(int i=1;i<=n;i++) {
for(int k=1;k<=n;k++) {
B.a[i][k]=B.a[i][k]*w[j];
if(i==k) B.a[i][k]++;
B.a[i][k]%=mod;
}
}
B=ksm(B,L);
f[j]=B.a[x][y];
}
for(int i=0;i<k;i++) F[i]=w[C_2(i)]*f[i]%mod;
for(int i=0;i<2*k;i++) G[i]=w[(k-C_2(i))%k]%mod;
reverse(G,G+2*k);
int d=ceil(log2(3*k));
mul(F,G,ans,d);
ll invk=ksm(k,mod-2);
for(int i=0;i<k;i++) {
cout<<ans[2*k-1-i]*invk%mod*w[C_2(i)]%mod<<"
";
}
return 0;
}