【笔记】有向图邻接矩阵幂敛指数及其周期的一种可行性算法
引入
一个 $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);
}