题目分析

构造邻接矩阵 $G$，要求的是 $\sum_{i=0}^n C_n^iG^i [k|i]$

单位根反演嘛，$\frac{1}{k}\sum_{i=0}^{k-1} \omega_k^{in}=[k|n]$。

$$\sum _ {i=0}^n C _ n^iG^i[k|i]=\frac{1}{k} \sum _ {i=0}^n C _ n^iG^i \sum _ {j=0}^{k-1} \omega^{ij} _ k=\frac{1}{k}\sum _ {j=0}^{k-1}\sum _ {i=0}^nC _ n^i\omega _ k^{ij}G^i=\frac{1}{k}\sum _ {j=0}^{k-1}(E+\omega _ k^jG)^n$$

据 boshi 说，找出 $P$的原根 $g$，单位根 $\omega_k=g^{\frac{P-1}{k}}$

找原根就是暴力找一个 $g$，若满足 $g^{P-1} \equiv 1 \pmod{P}$且对于 $P-1$的任何一个质因子 $p_i$，$g^{\frac{P-1}{p_i}} \not\equiv 1 \pmod{P}$，这个 $g$就是原根。

代码

#include<bits/stdc++.h>
using namespace std;
#define RI register int
typedef long long LL;
int K,P,m,w,S,T,road;LL n;
int pri[100005];
struct matrix{int t[5][5];}E,G,re,ans;

int ksm(int x,int y) {
    int re=1;
    for(;y;y>>=1,x=1LL*x*x%P) if(y&1) re=1LL*re*x%P;
    return re;
}
int qm(int x) {return x>=P?x-P:x;}
matrix operator * (matrix A,matrix B) {
    matrix C;
    for(RI i=0;i<m;++i)
        for(RI j=0;j<m;++j) C.t[i][j]=0;
    for(RI k=0;k<m;++k)
        for(RI i=0;i<m;++i)
            for(RI j=0;j<m;++j)
                C.t[i][j]=qm(C.t[i][j]+1LL*A.t[i][k]*B.t[k][j]%P);
    return C;
}
matrix operator * (int A,matrix B) {
    for(RI i=0;i<m;++i)
        for(RI j=0;j<m;++j) B.t[i][j]=1LL*B.t[i][j]*A%P;
    return B;
}
matrix operator + (matrix A,matrix B) {
    for(RI i=0;i<m;++i)
        for(RI j=0;j<m;++j) A.t[i][j]=qm(A.t[i][j]+B.t[i][j]);
    return A;
}
int getg(int x) {
    int kx=x-1,js=0;
    for(RI i=2;i*i<=kx;++i)
        if(kx%i==0) {pri[++js]=i;while(kx%i==0) kx/=i;}
    if(kx!=1) pri[++js]=kx;
    for(RI i=2;;++i) {
        int flag=1;
        if(ksm(i,x-1)!=1) continue;
        for(RI j=1;j<=js;++j) if(ksm(i,(x-1)/pri[j])==1) {flag=0;break;}
        if(flag) return i;
    }
}
matrix ksm_mt(matrix X,LL y) {
    matrix re=E;
    for(;y;y>>=1,X=X*X) if(y&1) re=re*X;
    return re;
}
int main()
{
    int x,y,ww;
    while(~scanf("%d%lld%d%d",&m,&n,&K,&P)) {
        w=ksm(getg(P),(P-1)/K),ww=1;
        for(RI i=0;i<m;++i)
            for(RI j=0;j<m;++j) E.t[i][j]=G.t[i][j]=ans.t[i][j]=0;
        for(RI i=0;i<m;++i) E.t[i][i]=1;
        scanf("%d%d%d",&road,&S,&T);
        for(RI i=1;i<=road;++i) scanf("%d%d",&x,&y),++G.t[x-1][y-1];
        for(RI i=0;i<K;++i) ans=ans+ksm_mt(E+ww*G,n),ww=1LL*ww*w%P;
        printf("%lld\n",1LL*ksm(K,P-2)*ans.t[S-1][T-1]%P);
    }
    return 0;
}