YbtOJ#723-欧拉之树【莫比乌斯反演,虚树】

1|0正题

题目链接:http://www.ybtoj.com.cn/contest/121/problem/2


1|1题目大意

给出n个点的一棵树,每个点有一个权值ai,求

i=1nj=1ndis(i,j)×φ(ai×aj)

2n2×105a恰好是一个排列。


1|2解题思路

一个十分显然的结论就是φ(x×y)=φ(x)×φ(y)×gcd(x,y)φ(gcd(x,y))。(相同的质因子只保留一个数p1的就好了)

然后顺便把点编号换一下使得ai=i再枚举约数就是

d=1nφ(d)di=1nj=1ndis(i,j)φ(i)φ(j)×[gcd(i,j)=d]

然后就可以莫反了,定义

gd=d|ind|jndis(i,j)φ(i)φ(j)

gd=d|ind|jn(depi+depj2deplca(i,j))φ(i)φ(j)

gd=2d|indepiφ(i)d|jnφ(j)2k=1ni=1nj=1n[lca(i,j)=k]depkφ(i)φ(j)

把所有d倍的点加入虚树,然后用树形dp计算后面那个东西,前面那个可以直接算。

然后答案就是

d=1nφ(d)dd|igi

时间复杂度O(nlog2n),有点卡常。


1|3code

#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__

本文作者QuantAsk
本文链接https://www.cnblogs.com/QuantAsk/p/14432453.html
关于博主:退役OIer,GD划水选手
版权声明:本博客所有文章除特别声明外,均采用 BY-NC-SA 许可协议。转载请注明出处!
声援博主:如果您觉得文章对您有帮助,可以点击文章右下角推荐一下。您的鼓励是博主的最大动力!
posted @   QuantAsk  阅读(51)  评论(0编辑  收藏  举报
编辑推荐:
· .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语句:使用策略模式优化代码结构
点击右上角即可分享
微信分享提示