题目链接:http://www.ybtoj.com.cn/contest/121/problem/2
给出n个点的一棵树,每个点有一个权值ai,求
n∑i=1n∑j=1dis(i,j)×φ(ai×aj)
2≤n≤2×105,a恰好是一个排列。
一个十分显然的结论就是φ(x×y)=φ(x)×φ(y)×gcd(x,y)φ(gcd(x,y))。(相同的质因子只保留一个数p−1的就好了)
然后顺便把点编号换一下使得ai=i再枚举约数就是
n∑d=1φ(d)dn∑i=1n∑j=1dis(i,j)φ(i)φ(j)×[gcd(i,j)=d]
然后就可以莫反了,定义
gd=n∑d|in∑d|jdis(i,j)φ(i)φ(j)
gd=n∑d|in∑d|j(depi+depj−2deplca(i,j))φ(i)φ(j)
gd=2n∑d|idepiφ(i)n∑d|jφ(j)−2n∑k=1n∑i=1n∑j=1[lca(i,j)=k]depkφ(i)φ(j)
把所有d倍的点加入虚树,然后用树形dp计算后面那个东西,前面那个可以直接算。
然后答案就是
n∑d=1φ(d)d∑d|igi
时间复杂度O(nlog2n),有点卡常。
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cctype>
#define ll long long
#pragma GCC optimize(2)
%:pragma GCC optimize(3)
%:pragma GCC optimize("Ofast")
%:pragma GCC optimize("inline")
using namespace std;
const int N=4e5+10,T=20,P=1e9+7;
int read() {
int 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-48,c=getchar();
return x*f;
}
struct node{
int to,next;
}a[N<<1];
int n,m,tot,top,dfc,ls[N],rfl[N];
int mu[N],phi[N],pri[N],prn,s[N],p[N];
int stn,lg[N],wz[N],rfn[N],dep[N],f[N][T];
ll S[N],dp[N],g[N],ans;bool v[N],mark[N];
ll power(ll x,ll b){
ll ans=1;
while(b){
if(b&1)ans=ans*x%P;
x=x*x%P;b>>=1;
}
return ans;
}
void prime(){
mu[1]=phi[1]=1;
for(int i=2;i<=n;i++){
if(!v[i])pri[++prn]=i,phi[i]=i-1,mu[i]=-1;
for(int j=1;j<=prn&&i*pri[j]<=n;j++){
v[i*pri[j]]=1;
if(i%pri[j]==0){
phi[i*pri[j]]=phi[i]*pri[j];
break;
}
phi[i*pri[j]]=phi[i]*(pri[j]-1);
mu[i*pri[j]]=-mu[i];
}
}
return;
}
void addl(int x,int y){
a[++tot].to=y;
a[tot].next=ls[x];
ls[x]=tot;return;
}
bool cmp(int x,int y)
{return rfn[x]<rfn[y];}
void dfs(int x,int fa){
dep[x]=dep[fa]+1;rfn[x]=++dfc;
f[++stn][0]=x;wz[x]=stn;
for(int i=ls[x];i;i=a[i].next){
int y=a[i].to;
if(y==fa)continue;
dfs(y,x);f[++stn][0]=x;
}
return;
}
int LCA(int l,int r){
l=wz[l];r=wz[r];
if(l>r)swap(l,r);
int z=lg[r-l+1],x=f[l][z],y=f[r-(1<<z)+1][z];
return dep[x]<dep[y]?x:y;
}
void Ins(int x){
if(!top){s[++top]=x;return;}
int lca=LCA(x,s[top]);
while(top>1&&dep[s[top-1]]>dep[lca])
addl(s[top-1],s[top]),top--;
if(dep[s[top]]>dep[lca])addl(lca,s[top]),top--;
if((!top)||s[top]!=lca)s[++top]=lca;
s[++top]=x;return;
}
void calc(int x,ll &ans){
if(mark[x])S[x]=phi[x],dp[x]=1ll*phi[x]*phi[x]%P;
else S[x]=dp[x]=0;
for(int i=ls[x];i;i=a[i].next){
int y=a[i].to;calc(y,ans);
(dp[x]+=S[x]*S[y]*2ll%P)%=P;
S[x]=(S[x]+S[y])%P;
}
(ans+=P-1ll*dp[x]*dep[x]%P)%=P;
ls[x]=mark[x]=0;return;
}
signed main()
{
freopen("sm.in","r",stdin);
freopen("sm.out","w",stdout);
n=read();prime();
for(int i=1;i<=n;i++){
int x=read();
rfl[i]=x;
}
for(int i=1;i<n;i++){
int x=read(),y=read();
x=rfl[x];y=rfl[y];
addl(x,y);addl(y,x);
}
dfs(1,1);
for(int j=1;(1<<j)<=stn;j++)
for(int i=1;i+(1<<j)-1<=stn;i++){
int x=f[i][j-1],y=f[i+(1<<j-1)][j-1];
f[i][j]=(dep[x]<dep[y])?x:y;
}
for(int i=2;i<=stn;i++)lg[i]=lg[i>>1]+1;
memset(ls,0,sizeof(ls));
for(int k=1;k<=n;k++){
m=top=tot=0;ll sum=0;
for(int i=k;i<=n;i+=k)
p[++m]=i,sum+=phi[i];
sort(p+1,p+1+m,cmp);sum%=P;
if(p[1]!=1)s[++top]=1;
for(int i=1;i<=m;i++){
Ins(p[i]);mark[p[i]]=1;
(g[k]+=1ll*phi[p[i]]*dep[p[i]]%P*sum%P)%=P;
}
while(top>1)addl(s[top-1],s[top]),top--;
calc(1,g[k]);g[k]=g[k]*2ll%P;
}
for(int i=1;i<=n;i++){
ll tmp=0;
for(int j=i;j<=n;j+=i)
(tmp+=mu[j/i]*g[j]%P)%=P;
(ans+=tmp*i%P*power(phi[i],P-2)%P)%=P;
}
printf("%d\n",(ans+P)%P*power(1ll*n*(n-1)%P,P-2)%P);
return 0;
}
__EOF__
【推荐】国内首个AI IDE,深度理解中文开发场景,立即下载体验Trae
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· .NET Core 中如何实现缓存的预热?
· 从 HTTP 原因短语缺失研究 HTTP/2 和 HTTP/3 的设计差异
· AI与.NET技术实操系列:向量存储与相似性搜索在 .NET 中的实现
· 基于Microsoft.Extensions.AI核心库实现RAG应用
· Linux系列:如何用heaptrack跟踪.NET程序的非托管内存泄露
· TypeScript + Deepseek 打造卜卦网站:技术与玄学的结合
· 阿里巴巴 QwQ-32B真的超越了 DeepSeek R-1吗?
· 【译】Visual Studio 中新的强大生产力特性
· 张高兴的大模型开发实战:(一)使用 Selenium 进行网页爬虫
· 【设计模式】告别冗长if-else语句:使用策略模式优化代码结构