Problem

OpenJ_POJ

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;
}