LOJ#3058. 「HNOI2019」白兔之舞 单位根反演+矩阵乘法+MTT
假如说每次只能跳一步，一共跳 $i$ 步到达 $(i,x)$ 的方案数是好求的：
$f[i][j]=\sum_{k=1}^{n} f[i-1][k] \times w(k,j)$.
时间复杂度为 $O(Ln^2)$.
$ans_{t}=\sum_{i=0}^{L} f_{L,y} \times \binom{L}{i} [i \mod k=t]$
考虑求解 $f_{L,y}$ 时运用矩阵乘法，然后对后面的特判进行单位根反演：
$ans_{t}=\sum_{i=0}^{L} F^i[x,y] \binom{L}{i}\frac{1}{k}\sum_{j=0}^{k-1}w_{k}^{j(i-t)}$
$ans_{t}=\frac{1}{k} \sum_{j=0}^{k-1} w_{k}^{-tj}\sum_{i=0}^{L}F^i[x,y]\binom{L}{i}w_{k}^{ij}$
后面部分可以写成二项式定理的形式，即 $(Fw_{k}^{j}+I)^L$
那么 $ans_{t}=\frac{1}{k} \sum_{j=0}^{k-1}w_{k}^{-tj}(Fw_{k}^{j}+I)^L$
令 $c_{j}=(Fw_{k}^{j}+I)^L[x,y]$，这个可以提前都预处理出来.
则 $ans_{t}=\frac{1}{k} \sum_{j=0}^{k-1} w_{k}^{-tj}c_{j}$
这里面 $-tj$ 是一个乘法的形式，没法再优化了.
有结论 $ij=\binom{i+j}{2}-\binom{i}{2}-\binom{j}{2}$
这样就没有相乘的形式了.
最后 $ans_{t}=\frac{w_{k}^{\binom{t}{2}}}{k}\sum_{j=0}^{k-1}w_{k}^{\binom{j}{2}}w_{k}^{-\binom{t+j}{2}}c_{j}$
这个可以直接用 MTT 来算.
但是这里面 $k$ 并不是 2 的整数次幂，但保证 $p-1$ 是 $k$ 的倍数，所以要求原根 $g$.
然后单位根就是 $g^{\frac{(p-1)}{k}}$.
根据费马小定理：$a^{p-1} \equiv 1(\mod p)$，$n$ 次单位根恰好为 $1$，符合要求.
求原根的话直接枚举原根，然后判断是否满足 $g^{\frac{p-1}{prime[i]}}$ 都不为 1 即可.
#include <vector>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
#define N 100009
#define ll long long
#define pb push_back
#define db long double
#define setIO(s) freopen(s".in","r",stdin)
using namespace std;
const db pi=acos(-1.0);
vector<int>g;
int gn,n,L,S,T,P,K,wi[N],val[5][5],seq[N<<1],os[N<<1];
int qpow(int x,int y,int mod) {
int tmp=1;
for(;y;y>>=1,x=(ll)x*x%mod) {
if(y&1) tmp=(ll)tmp*x%mod;
}
return tmp;
}
int getroot(int x) {
int phi=x-1;
for(int i=2;i*i<=x;++i) {
if(phi%i==0) {
g.pb(i);
while(phi%i==0) phi/=i;
}
}
if(phi>1) g.pb(phi);
phi=x-1;
for(int i=2;;++i) {
int flag=1;
for(int j=0;j<g.size()&&flag;++j) {
if(qpow(i,phi/g[j],x)==1) flag=0;
}
if(flag) return i;
}
}
ll sqr(ll x) {
return x*(x-1)/2;
}
void init() {
gn=qpow(getroot(P),(P-1)/K,P);
wi[0]=1;
for(int i=1;i<K;++i) {
wi[i]=(ll)wi[i-1]*gn%P;
}
}
struct Matrix {
int c[3][3];
Matrix(int x=0) {
memset(c,0,sizeof(c));
for(int i=0;i<3;++i) c[i][i]=x;
}
int *operator[](int x){ return c[x]; }
Matrix operator*(const Matrix &b) const {
Matrix an;
for(int i=0;i<3;++i)
for(int j=0;j<3;++j) {
for(int k=0;k<3;++k)
(an.c[i][j]+=(ll)c[i][k]*b.c[k][j]%P)%=P;
}
return an;
}
friend Matrix operator^(Matrix x,int y) {
Matrix an(1);
for(;y;y>>=1,x=x*x)
if(y&1) an=an*x;
return an;
}
}W;
const int base=1<<15;
struct Poly {
struct cp {
db x,y;
cp(db a=0,db b=0) { x=a,y=b; }
cp operator+(const cp &b) const { return cp(x+b.x,y+b.y); }
cp operator-(const cp &b) const { return cp(x-b.x,y-b.y); }
cp operator*(const cp &b) const { return cp(x*b.x-y*b.y,x*b.y+y*b.x); }
}f[2][N<<2],g[2][N<<2],ans[3][N<<2];
void FFT(cp *a,int len,int op) {
for(int i=0,k=0;i<len;++i) {
if(i>k) swap(a[i],a[k]);
for(int j=len>>1;(k^=j)<j;j>>=1);
}
for(int l=1;l<len;l<<=1) {
cp wn(cos(pi/l),op*sin(pi/l));
for(int i=0;i<len;i+=l<<1) {
cp w(1,0),x,y;
for(int j=0;j<l;++j) {
x=a[i+j],y=w*a[i+j+l];
a[i+j]=x+y,a[i+j+l]=x-y;
w=w*wn;
}
}
}
if(op==-1) {
for(int i=0;i<len;++i) a[i].x/=len;
}
}
ll actual(db x,int mod) {
return (ll)((ll)(x+0.5))%mod;
}
void MTT(int *a,int n,int *b,int m,int mod,int *c) {
int lim;
for(lim=1;lim<(n+m+1);lim<<=1);
for(int i=0;i<lim;++i) {
f[0][i]=f[1][i]=g[0][i]=g[1][i]=cp();
ans[0][i]=ans[1][i]=ans[2][i]=cp();
}
for(int i=0;i<=n;++i) {
f[0][i].x=a[i]/base;
f[1][i].x=a[i]%base;
}
for(int i=0;i<=m;++i) {
g[0][i].x=b[i]/base;
g[1][i].x=b[i]%base;
}
FFT(f[0],lim,1),FFT(f[1],lim,1);
FFT(g[0],lim,1),FFT(g[1],lim,1);
for(int i=0;i<lim;++i) {
ans[0][i]=f[0][i]*g[0][i];
ans[1][i]=f[0][i]*g[1][i]+f[1][i]*g[0][i];
ans[2][i]=f[1][i]*g[1][i];
}
for(int i=0;i<3;++i) {
FFT(ans[i],lim,-1);
}
for(int i=0;i<=n+m;++i) {
ll x=(ll)(actual(ans[0][i].x,mod)<<30ll)%mod;
ll y=(ll)(actual(ans[1][i].x,mod)<<15ll)%mod;
ll z=actual(ans[2][i].x,mod);
c[i]=(ll)((ll)(x+y)%mod+z)%mod;
}
}
}poly;
void get_matrix(int x) {
for(int i=0;i<3;++i) {
for(int j=0;j<3;++j) {
W[i][j]=(ll)val[i+1][j+1]*wi[x]%P;
}
++W[i][i];
}
}
int AA[N<<1],ANS[N*3];
int main() {
// setIO("input");
scanf("%d%d%d%d%d%d",&n,&K,&L,&S,&T,&P);
init();
for(int i=1;i<=n;++i) {
for(int j=1;j<=n;++j) scanf("%d",&val[i][j]);
}
for(int i=0;i<K;++i) {
get_matrix(i);
W=W^L;
seq[i]=(ll)W[S-1][T-1]*wi[(int)((ll)sqr(i)%K)]%P;
}
for(int i=0;i<=2*K-2;++i) {
os[i]=(ll)wi[(int)((ll)(K-(ll)sqr(i)%K+K)%K)];
}
for(int i=0;i<K;++i) AA[i]=seq[K-1-i];
poly.MTT(AA,K-1,os,2*K-2,P,ANS);
int inv=qpow(K,P-2,P);
for(int i=0;i<K;++i) {
printf("%d\n",(ll)inv*wi[(int)((ll)sqr(i)%K)]%P*ANS[K-1+i]%P);
}
return 0;
}
