Description

定义有向图邻接矩阵A的周期为最小的d，使得存在正整数k，对于任意n>=k，都有\(A^n=A^{n+d}\)

最小的k称为A的幂敛指数。

现给出一个n个点，m条边有向图，求它的邻接矩阵的周期对10^9+7取模的结果。

n<=100000,m<=200000

对于n<=200,m<=3000的数据，你还需要求出它的幂敛指数。

Solution

知乎上有一篇是讲这个玩意的，里面有一些参考文献（我也没看过），其中证明了一些结论。

k的上界是O(n^2)的

一个强联通图的d为其中所有的环长度的最大公约数

有向图的d为其中每个强联通分量的周期的最小公倍数

求d，我们可以先求出所有的强联通分量，对于每个强联通分量我们分别在上面dfs

若搜到了一条非树边（返祖边或横叉边）连接两个dfs树上深度为i,j的点，那么记d=|i-j+1|，要么存在一条

长为d的环，要么存在两个环长差为d

根据gcd(x,y)=gcd(x,x-y)，那么我们只需要将d与已经求得的答案取gcd即可

这样求一个有向图的周期的时间复杂度是\(O(m+n)\)的（排除了gcd和lcm的时间复杂度，这一部分可以分解质因数解决）

求k，k显然满足二分性，我们可以倍增找到最大的p，满足2^p<k，然后我们逐次将p-1，判断能否加上答案。

最后求出来的就是不满足周期性的最大的指数，+1就是k。

具体实现用矩阵乘法+快速幂，由于矩阵是01矩阵，因此可以bitset加速(枚举i,k，j可以整一行或过来)

时间复杂度\(O({n^3*\log d+n^3\log k\over \omega})\)

Code

#include 
    
     #define fo(i,a,b) for(int i=a;i<=b;++i)#define fod(i,a,b) for(int i=a;i>=b;--i)#define N 400005#define LL long longconst LL mo=1000000007;using namespace std;int n,m,tp,fs[N],nt[N],dt[N],m1,rs[N],ft[N],st[N],dfn[N],low[N],dep[N];void link(int x,int y){    nt[++m1]=fs[x];    dt[fs[x]=m1]=y;}int getf(int k){    return (ft[k]==k||ft[k]==0)?k:ft[k]=getf(ft[k]);}void merge(int x,int y){    int fx=getf(x),fy=getf(y);    if(fx!=fy) ft[fx]=fy;}int pr[N];bool bz[N];void prp(){    memset(bz,0,sizeof(bz));    fo(i,2,n)     {        if(!bz[i]) pr[++pr[0]]=i;        for(int j=1;j<=pr[0]&&i*pr[j]<=n;j++)         {            bz[i*pr[j]]=1;            if(i%j==0) break;        }    }}map
     
       hp;LL ksm(LL k,LL n){    LL s=1;    for(;n;n>>=1,k=k*k%mo) if(n&1) s=s*k%mo;    return s;}LL ans;void make(int k){    for(int i=1;pr[i]*pr[i]<=k;i++)    {        if(k%pr[i]==0)        {            int c=0;            while(k%pr[i]==0) c++,k/=pr[i];            if(!hp[pr[i]])             {                hp[pr[i]]=c;                if(!tp) ans=ans*ksm(pr[i],c)%mo;                else ans=ans*ksm(pr[i],c);            }            else            {                int v=hp[pr[i]];                if(c>v)                 {                    if(!tp) ans=ans*ksm(pr[i],c-v)%mo;                    else ans=ans*ksm(pr[i],c-v);                    hp[pr[i]]=c;                }            }        }    }    if(k>1&&!hp[k])     {        hp[k]=1;        if(!tp) ans=ans*k%mo;         else ans=ans*k;    }}void pop(int k){    while(st[st[0]]!=k)    {        bz[st[st[0]]]=0;        merge(st[st[0]],k);                     st[st[0]--]=0;    }    bz[k]=0,st[st[0]--]=0;}void tarjan(int k){    low[k]=dfn[k]=++dfn[0];    bz[st[++st[0]]=k]=1;    for(int i=fs[k];i;i=nt[i])     {        int p=dt[i];        if(!dfn[p]) tarjan(p),low[k]=min(low[k],low[p]);        else if(bz[p]) low[k]=min(low[k],dfn[p]);    }    if(low[k]>=dfn[k]) pop(k);}int gcd(int x,int y){       return(!y)?x:gcd(y,x%y);}void dfs(int k,int dp){    dep[k]=dp;    for(int i=fs[k];i;i=nt[i])    {        int p=dt[i];        if(getf(p)==getf(k))        {            if(!dep[p]) dfs(p,dp+1);            else             {                if(!rs[getf(k)]) rs[getf(k)]=abs(dep[k]+1-dep[p]);                else rs[getf(k)]=gcd(rs[getf(k)],abs(dep[k]+1-dep[p]));            }         }    }}struct mat{    bitset<201> a[201];    friend mat operator *(const mat &x,const mat &y)    {        mat z;        fo(i,1,n) z.a[i].reset();        fo(i,1,n) fo(k,1,n)             if(x.a[i][k]) z.a[i]|=y.a[k];        return z;    }    friend bool operator ==(const mat &x,const mat &y)    {        fo(i,1,n) if((x.a[i]^y.a[i]).any()) return 0;        return 1;     }}ap,wp[31];mat ksd(mat k,LL n){    n--;    mat s=k;    for(;n;n>>=1,k=k*k) if(n&1) s=s*k;    return s;}int main(){    cin>>n>>m>>tp;    fo(i,1,m)    {        int x,y;        scanf("%d%d",&x,&y);        link(x,y);        if(tp) ap.a[x][y]=1;    }    fo(i,1,n)         if(!dfn[i]) tarjan(i);    ans=1;    prp();    fo(i,1,n)     {        if(getf(i)==i)         {            rs[i]=0;            dfs(i,1);            make(rs[i]);        }    }    if(tp)    {        mat ws=ksd(ap,ans);        wp[0]=ap;        LL k=1,c=0;        while(!(wp[c]*ws==wp[c])) c++,k<<=1,wp[c]=wp[c-1]*wp[c-1];        k>>=1,c--;        mat sp=wp[c];        for(LL v=k>>1,c1=c-1;v;c1--,v>>=1)        {            mat np=sp*wp[c1];            if(!(np*ws==np)) sp=np,k+=v;         }        printf("%lld ",k+1);    }    printf("%lld\n",ans%mo);}