P8338-[AHOI2022]排列【质因数分解】

正题

题目链接:https://www.luogu.com.cn/problem/P8338

---

题目大意

给出一个排列\(p_i\)，定义\(a_i^0=i,a_i^k=a_{p_i}^{k-1}\)。

对于排列\(P\)，定义\(F(P)\)表示最小的一个正整数\(k\)满足\(P^{k+1}=P\)。
定义\(f(i,j)\)，若存在一个\(p_i^k=p_j\)那么\(f(i,j)=0\)，否则记\(P'\)表示将\(p_i\)和\(p_j\)交换后的排列，\(f_{i,j}=F(P')\)

求

\[\sum_{i=1}^n\sum_{j=1}^nf(i,j) \]

答案对\(10^9+7\)取模

\(1\leq n\leq 5\times 10^5,1\leq T\leq 5\)

---

解题思路

考虑置换环，\(i\rightarrow p_i\)，那么\(F(P)\)就是所有环的大小的LCM。
然后交换\(p_i,p_j\)的话就相当于把两个环拼到一起。

注意到所有环的长度和为\(n\)，那么就证明不同的长度不超过\(\sqrt n\)种，我们可以考虑同一种长度一起处理。

因为最多删去两个环，处理出每个环长质因数幂数最大的前三个。

然后暴力处理就好了，需要记得提前预处理处所有数的质因数分解。

时间复杂度：\(O(Tn\log n)\)

---

code

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<vector>
#include<cctype>
#include<stack>
#define ll long long
#define mp(x,y) make_pair(x,y)
using namespace std;
const ll N=5e5+10,P=1e9+7;
ll T,n,m,ans,p[N],c[N];
ll inv[N],val[N],fa[N];
bool v[N],vis[N];
pair<ll,ll> a[N];
vector<pair<ll,ll> >q[N];
vector<ll>z[N],pw[N];
stack<ll> cl;
ll read(){
	ll x=0,f=1;char c=getchar();
	while(!isdigit(c)){if(c=='-')f=-f;c=getchar();}
	while(isdigit(c)){x=(x<<1)+(x<<3)+c-'0';c=getchar();}
	return x*f;
}
ll dfs(ll x){
	if(v[x])return 0;v[x]=1;
	return dfs(fa[x])+1;
}
void Ins(ll x){
	for(ll i=0;i<q[x].size();i++){
		ll a=q[x][i].first,b=q[x][i].second;
		z[a].push_back(b);
	}
}
void Del(ll x){
	for(ll i=0;i<q[x].size();i++){
		ll a=q[x][i].first,b=q[x][i].second;
		if(b==z[a][p[a]]){
			ans=ans*inv[pw[a][z[a][p[a]]]]%P;
			if(!p[a])cl.push(a);p[a]++;
			ans=ans*pw[a][z[a][p[a]]]%P;
		}
	}
	return;
}
void Add(ll x){
	for(ll i=0;i<q[x].size();i++){
		ll a=q[x][i].first,b=q[x][i].second;
		if(b>z[a][p[a]]){
			ans=ans*inv[pw[a][z[a][p[a]]]]%P;
			ans=ans*pw[a][b]%P;
		}
	}
}
bool cmp(ll x,ll y){return x>y;}
signed main()
{
//	freopen("perm3.in","r",stdin);
	T=read();inv[0]=inv[1]=1;
	for(ll i=2;i<N;i++)inv[i]=P-inv[P%i]*(P/i)%P;
	for(ll i=1;i<N;i++)val[i]=i;
	for(ll i=2;i<N;i++){
		if(!vis[i]){
			pw[i].push_back(1);
			for(ll j=i;j<N;j=j*i)
				pw[i].push_back(j);
			for(ll j=i+i;j<N;j+=i)vis[j]=1;
		}
		for(ll j=i;j<N;j+=i){
			if(val[j]%i==0){
				ll c=0;
				while(val[j]%i==0)val[j]/=i,c++;
				q[j].push_back(mp(i,c));
			}
		}
	}
	while(T--){
		m=0;n=read();
		for(ll i=1;i<=n;i++)
			z[i].clear(),p[i]=v[i]=c[i]=0;
		for(ll i=1;i<=n;i++)fa[i]=read();
		for(ll i=1;i<=n;i++)
			if(!v[i])c[dfs(i)]++;
		for(ll i=1;i<=n;i++)
			if(c[i]){
				a[++m]=mp(i,c[i]);
				for(ll j=1;j<=c[i];j++)Ins(i);
			}
		ans=1;
		for(ll i=2;i<=n;i++)
			if(!vis[i]){
				z[i].push_back(0);
				sort(z[i].begin(),z[i].end(),cmp);
				p[i]=0;ans=ans*pw[i][z[i][0]]%P;
			}
		ll pre=ans,sum=0;
		for(ll i=1;i<=m;i++){
			for(ll j=1;j<i;j++){
				Del(a[i].first);Del(a[j].first);
				Add(a[i].first+a[j].first);
				(sum+=2ll*ans*a[i].second*a[i].first%P*a[j].second*a[j].first%P)%=P;
				while(!cl.empty())p[cl.top()]=0,cl.pop();ans=pre;
			}
			if(a[i].second>1){
				Del(a[i].first);Del(a[i].first);
				Add(a[i].first*2);
				(sum+=ans%P*a[i].second*a[i].first%P*(a[i].second-1)*a[i].first%P)%=P;
				while(!cl.empty())p[cl.top()]=0,cl.pop();ans=pre;
			}
		}
		printf("%lld\n",sum);
	}
	return 0;
}