Problem
Solution
如果我们用 \(G\) 来表示邻接矩阵，那么答案其实就是求\(\sum_{k|i}^n \binom n i G^i\)
为了消除整除的限制，我们可以考虑单位根反演：
\[\frac 1 k\sum_{i=0}^{k-1}w_k^{in}=[k|n]
\]
\]
那么原题就是要求
\[\frac 1 k\sum_{i=0}^n \binom n i G^i\sum_{j=0}^{k-1}w_k^{ij}
\]
\]
\[\frac 1 k\sum_{j=0}^{k-1}\sum_{i=0}^n \binom n i (Gw_k^j)^i
\]
\]
\[\frac 1 k\sum_{j=0}^{k-1}(Gw_k^j+1)^n
\]
\]
注意到题目的一个限制\(P\equiv 1 \pmod k\)，此时有\(w_k=g^{\frac {p-1} {k}}\)。为什么呢？因为\(w_k^k=1=g^{p-1}\)，而且显然根据原根的性质，\(\forall_{i<k} w_k^i\)都两两不相等。
时间复杂度\(O(km^3\log n)\)
Code
#include <cstring>
#include <cstdio>
using namespace std;
typedef long long ll;
template <typename Tp> inline int getmin(Tp &x,Tp y){return y<x?x=y,1:0;}
template <typename Tp> inline int getmax(Tp &x,Tp y){return y>x?x=y,1:0;}
template <typename Tp> inline void read(Tp &x)
{
x=0;int f=0;char ch=getchar();
while(ch!='-'&&(ch<'0'||ch>'9')) ch=getchar();
if(ch=='-') f=1,ch=getchar();
while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=getchar();
if(f) x=-x;
}
int m,k,l,s,t,cnt,mod,G,wk,ans,pri[40];
ll n;
int pls(int x,int y){return x+y>=mod?x+y-mod:x+y;}
int dec(int x,int y){return x-y<0?x-y+mod:x-y;}
struct Matrix{
int a[5][5];
Matrix(){memset(a,0,sizeof(a));}
int* operator [] (int x){return a[x];}
void clear(){memset(a,0,sizeof(a));}
Matrix operator * (Matrix &b)const
{
Matrix res;
for(int i=0;i<5;i++)
for(int k=0;k<5;k++)
for(int j=0;j<5;j++)
res[i][j]=pls(res[i][j],(ll)a[i][k]*b[k][j]%mod);
return res;
}
Matrix operator + (Matrix b)const
{
Matrix res;
for(int i=0;i<5;i++)
for(int j=0;j<5;j++)
res[i][j]=pls(a[i][j],b[i][j]);
return res;
}
Matrix operator * (const int &b)const
{
Matrix res;
for(int i=0;i<5;i++)
for(int j=0;j<5;j++)
res[i][j]=(ll)a[i][j]*b%mod;
return res;
}
}A,B,C,E;
int power(int x,int y)
{
int res=1;
for(;y;y>>=1,x=(ll)x*x%mod)
if(y&1ll)
res=(ll)res*x%mod;
return res;
}
Matrix mat_pow(Matrix x,ll y)
{
Matrix res=x;
for(--y;y;y>>=1,x=x*x)
if(y&1)
res=res*x;
return res;
}
void find_G()
{
int x=mod-1,f;
for(int i=2;i*i<=x;i++)
if(x%i==0)
{
pri[++cnt]=i;
while(x%i==0) x/=i;
}
if(x>1) pri[++cnt]=x;
for(G=2;;G++)
{
f=1;
for(int j=1;j<=cnt;j++)
if(power(G,(mod-1)/pri[j])==1)
{f=0;break;}
if(f) return ;
}
}
int input()
{
int u,v;
if(scanf("%d",&m)==EOF) return 0;
cnt=0;A.clear();B.clear();C.clear();
read(n);read(k);read(mod);
read(l);read(s);read(t);--s;--t;
for(int i=1;i<=l;i++)
{
read(u);read(v);
++C[v-1][u-1];
}
for(int i=0;i<5;i++) E[i][i]=1;
find_G();B[s][0]=1;
return 1;
}
int main()
{
while(input())
{
wk=power(G,(mod-1)/k);
for(int i=0;i<k;i++,C=C*wk)
A=A+mat_pow(C+E,n);
A=A*B;
ans=(ll)A[t][0]*power(k,mod-2)%mod;
printf("%d\n",ans);
}
return 0;
}