题目概述

LOJ3058

解题报告

这个 $m\bmod k=t$ 以及 $k$ 是 $MOD-1$ 因子显然是单位根反演，所以考虑用生成函数表示答案。

令转移矩阵为 $M$ ，定义 $F_i(x)$ 表示第一维前 $i$ 个的生成函数，$x^j$ 的系数为走了 $j$ 步得到的所有矩阵之和，则：

$$ F_i(x)=F_{i-1}(x)+F_{i-1}(x)Mx=F_{i-1}(x)(Mx+I)=\cdots=(Mx+I)^{i} $$

最终生成函数为 $F(x)=(Mx+I)^{L}$ ，接下来套用单位根反演：

$$ ans_t={1\over k}\sum_{i=0}^{k-1}{F(\omega_k^i)\over\omega_k^{it}}={1\over k}\sum_{i=0}^{k-1}\omega_k^{-it}(M\omega_k^i+I)_{x,y}^{L} $$

由于要求所有的 $ans_t$ ，所以需要考虑卷积，把 $\omega_k^{-it}$ 用Bluestein套路换掉：$-it=-{i+t\choose 2}+{i\choose 2}+{t\choose 2}$ 。

$$ ans_t={\omega_k^{t\choose 2}\over k}\sum_{i=0}^{k-1}\omega_k^{-{i+t\choose 2}}\cdot \omega_k^{i\choose 2}(M\omega_k^i+I)_{x,y}^{L} $$

这样就是卷积形式了，套个MTT就行。

示例程序

#include<cmath>
#include<cstdio>
#include<algorithm>
#define fr first
#define sc second
#define mp make_pair
using namespace std;
typedef long long LL;typedef double DB;
typedef pair<DB,DB> CN;
const int maxs=32767,maxk=65536,maxt=1<<18;
const DB pi=3.1415926535897932384626433832795;

int n,K,L,S,T,MOD,gn,ig;
int p[maxs+5];bool pri[maxs+5];
struct Matrix{
    int s[3][3];
    void zero() {for (int i=0;i<n;i++) for (int j=0;j<n;j++) s[i][j]=0;}
    void unit() {zero();for (int i=0;i<n;i++) s[i][i]=1;}
}M,tem;
CN wn[maxt+5];int F[maxt+5],G[maxt+5];

__attribute__((constructor)) void Pre(){
    for (int i=0;i<(maxt>>1);i++) wn[i+(maxt>>1)]=mp(cos(2*pi*i/maxt),sin(2*pi*i/maxt));
    for (int i=(maxt>>1)-1;i;i--) wn[i]=wn[i<<1];
    for (int i=2;i<=maxs;i++){
        if (!pri[i]) p[++p[0]]=i;
        for (int j=1,t;j<=p[0] && (t=i*p[j])<=maxs;j++)
            {pri[t]=true;if (!(i%p[j])) break;}
    }
}
inline int ADD(int x,int y) {return x+y>=MOD?x+y-MOD:x+y;}
inline int MUL(int x,int y) {return (LL)x*y%MOD;}
int Pow(int w,LL b) {int s;for (s=1;b;b>>=1,w=MUL(w,w)) if (b&1) s=MUL(s,w);return s;}
int getg(int MOD){
    for (int i=2;i<MOD;i++){
        for (int j=1;j<=p[0] && (LL)p[j]*p[j]<=MOD;j++)
            if (!((MOD-1)%p[j])) if (Pow(i,(MOD-1)/p[j])==1) goto GG;
        return i;GG:continue;
    }
    return -1;
}
Matrix operator * (const Matrix &a,const Matrix &b){
    static Matrix c;c.zero();
    for (int i=0;i<n;i++)
        for (int j=0;j<n;j++)
            for (int k=0;k<n;k++)
                c.s[i][j]=ADD(c.s[i][j],MUL(a.s[i][k],b.s[k][j]));
    return c;
}
Matrix Pow(Matrix w,int b){
    static Matrix s;
    for (s.unit();b;b>>=1,w=w*w) if (b&1) s=s*w;
    return s;
}
inline CN operator ! (const CN &a) {return mp(a.fr,-a.sc);}
inline CN operator + (const CN &a,const CN &b) {return mp(a.fr+b.fr,a.sc+b.sc);}
inline CN operator - (const CN &a,const CN &b) {return mp(a.fr-b.fr,a.sc-b.sc);}
inline CN operator * (const CN &a,const CN &b) {return mp(a.fr*b.fr-a.sc*b.sc,a.fr*b.sc+a.sc*b.fr);}
void FFT(CN *a,int n,int f){
    CN x,y;
    if (f>0){
        for (int k=n>>1;k;k>>=1)
            for (int i=0;i<n;i+=k<<1)
                for (int j=0;j<k;j++){
                    x=a[i+j];y=a[i+j+k];
                    a[i+j+k]=(x-y)*wn[k+j];
                    a[i+j]=x+y;
                }
    } else {
        for (int k=1;k<n;k<<=1)
            for (int i=0;i<n;i+=k<<1)
                for (int j=0;j<k;j++){
                    x=a[i+j];y=a[i+j+k]*!wn[k+j];
                    a[i+j+k]=x-y;
                    a[i+j]=x+y;
                }
        for (int i=0;i<n;i++) a[i].fr/=n,a[i].sc/=n;
    }
}
void MTT(int *c,int *a,int n,int *b,int m){
    static int lg=15,BA=(1<<lg)-1;
    static CN A[maxt+5],P[maxt+5],Q[maxt+5];
    int t;for (t=1;t<=n+m;t<<=1);
    for (int i=0;i<=n;i++) P[i]=mp(a[i]&BA,a[i]>>lg);
    for (int i=n+1;i<t;i++) P[i]=mp(0,0);
    for (int i=0;i<=m;i++) Q[i]=mp(b[i]&BA,b[i]>>lg);
    for (int i=m+1;i<t;i++) Q[i]=mp(0,0);
    FFT(P,t,1);FFT(Q,t,1);
    for (int i=0;i<t;i++) A[i]=P[i];
    P[0]=mp((A[0].fr+A[0].fr)/2,(A[0].sc-A[0].sc)/2)*Q[0];
    Q[0]=mp((A[0].sc+A[0].sc)/2,(A[0].fr-A[0].fr)/2)*Q[0];
    for (int k=1,i=1,j;k<t;k<<=1)
        for (j=i^k-1;i<(k<<1);i++,j=i^k-1)
            P[i]=mp((A[i].fr+A[j].fr)/2,(A[i].sc-A[j].sc)/2)*Q[i],
            Q[i]=mp((A[i].sc+A[j].sc)/2,(A[j].fr-A[i].fr)/2)*Q[i];
    FFT(P,t,-1);FFT(Q,t,-1);
    for (int i=0;i<=n+m;i++){
        LL A0B0=LL(P[i].fr+0.5)%MOD,A0B1=LL(P[i].sc+0.5)%MOD;
        LL A1B0=LL(Q[i].fr+0.5)%MOD,A1B1=LL(Q[i].sc+0.5)%MOD;
        c[i]=(A0B0+(A0B1+A1B0<<lg)%MOD+(A1B1<<lg+lg)%MOD)%MOD;
    }
}
int main(){
    scanf("%d%d%d%d%d%d",&n,&K,&L,&S,&T,&MOD);S--;T--;
    gn=Pow(getg(MOD),(MOD-1)/K);ig=Pow(gn,MOD-2);
    for (int i=0;i<n;i++)
        for (int j=0;j<n;j++)
            scanf("%d",&M.s[i][j]);
    for (int i=0;i<=(K-1<<1);i++) F[i]=Pow(ig,(LL)i*(i-1)>>1);
    for (int t=0,pw=1;t<K;t++,pw=MUL(pw,gn)){
        tem.unit();
        for (int i=0;i<n;i++)
            for (int j=0;j<n;j++)
                tem.s[i][j]=ADD(tem.s[i][j],MUL(M.s[i][j],pw));
        tem=Pow(tem,L);
        G[K-1-t]=MUL(tem.s[S][T],Pow(gn,(LL)t*(t-1)>>1));
    }
    MTT(F,F,K-1<<1,G,K-1);
    for (int i=0,INV=Pow(K,MOD-2);i<K;i++){
        int pw=Pow(gn,(LL)i*(i-1)>>1);
        printf("%d\n",MUL(F[i+K-1],MUL(pw,INV)));
    }
    return 0;
}