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^3log kover 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
",ans%mo);
}