LOJ#3058. 「HNOI2019」白兔之舞 单位根反演+矩阵乘法+MTT

假如说每次只能跳一步，一共跳 $i$ 步到达 $(i,x)$ 的方案数是好求的：

$f[i][j]=sum_{k=1}^{n} f[i-1][k] imes w(k,j)$.

时间复杂度为 $O(Ln^2)$.

$ans_{t}=sum_{i=0}^{L} f_{L,y} imes inom{L}{i} [i mod k=t]$

考虑求解 $f_{L,y}$ 时运用矩阵乘法，然后对后面的特判进行单位根反演：

$ans_{t}=sum_{i=0}^{L} F^i[x,y] inom{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]inom{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=inom{i+j}{2}-inom{i}{2}-inom{j}{2}$

这样就没有相乘的形式了.

最后 $ans_{t}=frac{w_{k}^{inom{t}{2}}}{k}sum_{j=0}^{k-1}w_{k}^{inom{j}{2}}w_{k}^{-inom{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
",(ll)inv*wi[(int)((ll)sqr(i)%K)]%P*ANS[K-1+i]%P);   
    }
    return 0;
}