【笔记】有向图邻接矩阵幂敛指数及其周期的一种可行性算法

【笔记】有向图邻接矩阵幂敛指数及其周期的一种可行性算法

引入

一个 $n$ 个点有向图的邻接矩阵 $A$ 是一个 $n$ 阶布尔矩阵,它的幂组成的序列具有周期性。

设 $k,d$ 是满足 $A^k=A^{k+d}$ 的最小正整数,则 $k$ 称为 $A$ 的幂敛指数(index),$d$ 称为 $A$ 的周期(period)。

对于如下邻接矩阵 $A$,其各次幂分别为

$$A^1 = \left[\begin{matrix}0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 & 0 \end{matrix}\right]$$

$$A^2 = \left[ \begin{array}{ccccc} 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \end{array} \right]$$

$$A^3 = \left[ \begin{array}{ccccc} 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \end{array} \right]$$

$$A^4 = \left[ \begin{array}{ccccc} 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 & 0 \end{array} \right]$$

$$A^5 = \left[ \begin{array}{ccccc} 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \end{array} \right]$$

此时,$A^5 = A^2$,故 $k=2,d=3$。

思路

矩阵 $A$ 的周期等于在图上每个点经过一定步数所能到达的点的周期。根据 B. De Schutter and B. De Moor 的报告 On the sequence of consecutive powers of a matrix in a Boolean algebra,强连通图的周期是其包含的所有环长度的最大公约数,而一般图的周期是其所有强连通分量周期的最小公倍数。简单证明如下:

对于一个强连通分量,设其包含的所有环长度的最大公约数为 $d_0$,那么任意选取一点为起点,沿有向边循坏按 $1-d_0$ 标号,每个点都只会有唯一的标号。那么,经过无限步后,一个点能够到达的点集合一定是以 $d_0$ 为周期循环的。

那么,整张图的周期,显然就是每个强连通分量周期的最小公倍数。得到周期 $d$ 之后,幂敛指数 $k$ 可以较为轻松地计算求得。

算法实现

计算强连通分量使用 $\text{Tarjan}$ 算法可以获得 $\Theta(V+E)$ 的复杂度。

对于每个强连通分量,任取一个点作为根开始 $\text{DFS}$。过程中,一条两端深度分别为 $i,j$ 的非树边对答案的贡献 $d’=|i+1-j|$,代表着一个长度为 $d’$ 的环,或者一对长度差为 $d’$ 的环。找到所有的 $d’$,它们的最大公约数就是这个强连通分量的周期。

得到整张图的周期后,从 $p=0$ 开始,我们不断验证 $A^{2^p} = A^{2^p+d}$ 这一等式的合法性,每次将 $p+1$。$k$ 的范围为 $[2^{p-1},2^p]$,其中 $p$ 为等式第一次合法时的值。在这个区间内,直接二分查找使 $A^k=A^{k+d}$ 存在的最小值。因为 $d$ 比较大,确定 $k$ 的范围的第一次查找需要计算 $\Theta(\log d)$ 次矩阵乘法。使用 $\text{bitset}$ 优化矩阵乘法可以获得 $\Theta(n^3\log d/32)$ 的复杂度。

代码

#include <bits/stdc++.h>
using namespace std;
const int maxn=1e5+1;
const int maxl=2e2+1;
const int md=1e9+7;
typedef long long ll;

int read(){
	int x=0;char ch=getchar();
	while(!isdigit(ch)) ch=getchar();
	while(isdigit(ch)) x=(x<<3)+(x<<1)+ch-'0',ch=getchar();
	return x;
}

ll gcd(ll x,ll y){
	return y?gcd(y,x%y):x;
}

ll lcm(ll x,ll y){
	return y?x*y/gcd(x,y):x;
}

vector<int>to[maxn];
int n,m,ty;

int dfn[maxn],low[maxn],id[maxn];
int vis[maxn],sta[maxn];
ll nd[maxn],ans_d=1;
int val[maxn];
int top,cnt,tot;

void Tarjan(int u){
	dfn[u]=low[u]=++cnt;
	sta[++top]=u;vis[u]=1;
	for (auto v:to[u]){
		if (!dfn[v]){
			Tarjan(v);
			low[u]=min(low[u],low[v]);
		}
		else if (vis[v]) low[u]=min(low[u],dfn[v]);
	}
	if (dfn[u]==low[u]){
		++tot;
		int k=-1;
		while(k!=u){
			k=sta[top--];
			id[k]=tot;
			vis[k]=0;
		}
	}
}

void dfs(int u,int d){
    vis[u]=1;val[u]=d;
    for (auto v:to[u]){
        if (id[v]!=id[u]) continue; 
        if (!vis[v]) dfs(v,d+1);
		else if (val[v]!=d+1){
            nd[id[u]]=gcd(nd[id[u]],(ll)abs(d+1-val[v]));
        }
    }
}

struct matrix{
	bitset<maxl>a[maxl];
	inline void set(){
		for (register int i=1;i<=n;++i) a[i][i]=1;
	}
	inline matrix operator * (const matrix x) const{
		matrix res;
		for (register int i=1;i<=n;++i)
		for (register int k=1;k<=n;++k) if (a[i][k])
		res.a[i]|=x.a[k];
		return res;
	}
	inline bool operator == (const matrix x) const{
		for (register int i=1;i<=n;++i)
		if (a[i]!=x.a[i]) return 0;
		return 1;
	}
}grp,sone;

matrix qpow(matrix a,ll b){
	matrix res=sone;
	while(b){
		if (b&1) res=res*a;
		b>>=1;a=a*a;
	}
	return res;
}

ll qpow(ll a,ll b){
	ll res=1;
	while(b){
		if (b&1) res=res*a%md;
		a=a*a%md;b>>=1;
	}
	return res;
}

ll solve(){
	int p=1;
	matrix lst=grp;
	while(1){
		if (qpow(grp,p+ans_d)==lst) break;
		lst=qpow(lst,2);p<<=1;
	}
	ll l=max(p/2,1),r=p,res=1;
	while(l<=r){
		ll mid=(l+r)>>1;
		if (qpow(grp,mid+ans_d)==qpow(grp,mid))
		res=mid,r=mid-1;
		else l=mid+1;
	}
	return res;
}

int prm[maxn];

int main(){
	n=read(),m=read();
	sone.set();
	for (register int i=1;i<=m;++i){
		int a=read(),b=read();
		to[a].push_back(b);
		grp.a[a][b]=1;
	}
	for (register int i=1;i<=n;++i)
	if (!dfn[i]) Tarjan(i);
	memset(vis,0,sizeof(vis));
	for (register int i=1;i<=n;++i)
	if (!vis[i]) dfs(i,0);
	for (register int i=1;i<=tot;++i){
		for (register int k=2;k<=nd[i];++k){
			if (nd[i]%k) continue;
			int tmp=0;
			while(!(nd[i]%k)) tmp++,nd[i]/=k;
			prm[k]=max(prm[k],tmp);
		}
	}
	for (register int i=1;i<=n;++i)
		ans_d*=qpow(i,prm[i]);
	printf("%lld ",solve());
	printf("%lld\n",ans_d);
}