有向图邻接矩阵的幂敛指数与周期【图论】

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 <bits/stdc++.h>
#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 long
const 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<int,int> 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);
}

转载于:https://www.cnblogs.com/BAJimH/p/11018158.html