zoukankan      html  css  js  c++  java
  • Codeforces.809E.Surprise me!(莫比乌斯反演 虚树)

    题目链接

    (Description)

      给定一棵树,求$$frac{1}{n(n-1)/2} imessum_{iin[1,n],jin[1,n],i eq j}varphi(a_i imes a_j) imes dis(i,j) (mod 10^9+7)$$
      其中(a_i)([1,n])的一个排列,两两不同。

    (Solution)

      前面直接最后乘逆元就可以。看后面的(sum)怎么化。
      要想办法把(varphi(a_i imes a_j))化掉。
      考虑(varphi)的定义式:(varphi(n)=n imesprod_{i=1}^k(1-frac{1}{p_i}))
      (varphi(a_i)=a_i imesprod_{l=1}^k(1-frac{1}{p_l}))
      (varphi(a_j)=a_j imesprod_{l=1}^{k'}(1-frac{1}{p'_l}))
      两者相乘的话会多算(g=gcd(a_i,a_j))(prod_{l=1}^{k''}(1-frac{1}{p''_l}))部分。所以有:$$varphi(a imes b)=varphi(a) imesvarphi(b) imesfrac{gcd(a,b)}{varphi(gcd(a,b))}$$
      (sum)就可以化成上面那种形式。
      带着(gcd(a,b)),想到去枚举(d=gcd(a,b)),于是$$原式=sum_{d=1}^nfrac{d}{varphi(d)}sum_{i=1}^nsum_{j=1}^nvarphi(a_i)varphi(a_j)dis(i,j)[gcd(a_i,a_j)=d]$$
      令(f(d)=sum_{i=1}^nsum_{j=1}^nvarphi(a_i)varphi(a_j)dis(i,j)[gcd(a_i,a_j)=d])
      (F(d)=sum_{i=1}^nsum_{j=1}^nvarphi(a_i)varphi(a_j)dis(i,j)[dmid gcd(a_i,a_j)])
      那么要求的$$f(d)=sum_{dmid n}mu(frac{n}{d})F(n)=sum_{i=1}^{lfloorfrac{n}{d} floor}mu(i)F(id)$$
      如果能求出(F(id)),那么答案就可以在(O(nln n))时间内得到了。

      (F(d))中要求(dmidgcd(a_i,a_j)),那么只有(dmid a_i)的点之间才会产生贡献,于是我们可以建虚树。(d)最大是(n),总点数是(nln n)的。
      现在考虑如何在虚树上求答案。(dis(i,j))可以拆掉,写成:(dep(i)+dep(j)-2*dep(lca(i,j)))
      现在$$F(d)=sum_{1leq ileq n,dmid a_i}sum_{1leq jleq n,dmid a_j}varphi(a_i)varphi(a_j)[dep(i)+dep(j)-dep(lca(i,j))]$$
      (dep(i))(dep(j))可以在处理(i,j)时分别直接算出来,现在要算$$-2 imessum_{1leq ileq n,dmid a_i}sum_{1leq jleq n,dmid a_j}varphi(a_i)varphi(a_j)dep(lca(i,j))$$
      枚举lca,求其子树(varphi)的和就行了。
      构建虚树时 按照DFS序存下每个点,这样就可以省掉排序了。


      改了一个小时的MLE,发现:为什么数组开小MLE了==?而且N至少是3e5不能是2e5why。。


    //1107ms	75700KB
    #include <cmath>
    #include <cstdio>
    #include <cctype>
    #include <vector>
    #include <cstring>
    #include <algorithm>
    //#define gc() getchar()
    #define MAXIN 300000
    #define gc() (SS==TT&&(TT=(SS=IN)+fread(IN,1,MAXIN,stdin),SS==TT)?EOF:*SS++)
    #define pb push_back
    #define Vec std::vector<int>
    #define mod (1000000007)
    #define Add(x,v) (x)+=(v),(x)>=mod&&((x)-=mod)
    typedef long long LL;
    const int N=300005;//2e5+5???
    
    int n,A[N],Enum,H[N],nxt[N<<1],to[N<<1],dfn[N],Index,fa[N],sz[N],dep[N],son[N],tp[N],sk[N],top;
    int cnt,P[N],phi[N],mu[N],sum[N],f[N],inv[N];
    bool not_P[N];
    Vec fac[N],pt[N];
    char IN[MAXIN],*SS=IN,*TT=IN;
    
    inline int read()
    {
    	int now=0;register char c=gc();
    	for(;!isdigit(c);c=gc());
    	for(;isdigit(c);now=now*10+c-'0',c=gc());
    	return now;
    }
    inline void Add_direct(int u,int v){
    	to[++Enum]=v, nxt[Enum]=H[u], H[u]=Enum;
    }
    inline void AddEdge(int u,int v)
    {
    	to[++Enum]=v, nxt[Enum]=H[u], H[u]=Enum;
    	to[++Enum]=u, nxt[Enum]=H[v], H[v]=Enum;
    }
    inline int LCA(int u,int v)
    {
    	while(tp[u]!=tp[v]) dep[tp[u]]>dep[tp[v]]?u=fa[tp[u]]:v=fa[tp[v]];
    	return dep[u]>dep[v]?v:u;
    }
    void Pre()
    {
    	mu[1]=phi[1]=1;
    	for(int i=2; i<=n; ++i)
    	{
    		if(!not_P[i]) P[++cnt]=i, mu[i]=-1, phi[i]=i-1;
    		for(int j=1,v; j<=cnt&&(v=i*P[j])<=n; ++j)
    		{
    			not_P[v]=1;
    			if(i%P[j]) mu[v]=-mu[i], phi[v]=phi[i]*(P[j]-1);
    			else {mu[v]=0, phi[v]=phi[i]*P[j]; break;}
    		}
    	}
    	for(int i=1; i<=n; ++i)
    		for(int j=i; j<=n; j+=i) fac[j].push_back(i);
    	inv[1]=1;
    	for(int i=2; i<=n; ++i) inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
    }
    void preDFS1(int x)
    {
    	const Vec &fc=fac[A[x]];
    	for(int i=0,l=fc.size(); i<l; ++i)  pt[fc[i]].push_back(x);
    	int mx=0; sz[x]=1, dfn[x]=++Index;
    	for(int v,i=H[x]; i; i=nxt[i])
    		if((v=to[i])!=fa[x])
    		{
    			fa[v]=x, dep[v]=dep[x]+1, preDFS1(v), sz[x]+=sz[v];
    			if(sz[v]>mx) mx=sz[v], son[x]=v;
    		}
    }
    void preDFS2(int x,int _tp)
    {
    	tp[x]=_tp;
    	if(son[x])
    	{
    		preDFS2(son[x],_tp);
    		for(int i=H[x]; i; i=nxt[i])
    			if(to[i]!=son[x]&&to[i]!=fa[x]) preDFS2(to[i],to[i]);
    	}
    }
    void Insert(int x)
    {
    	if(top==1) {sk[++top]=x; return;}
    	int lca=LCA(x,sk[top]);
    	while(dfn[sk[top-1]]>=dfn[lca]) Add_direct(sk[top],sk[top--]);
    	if(sk[top]!=lca) Add_direct(lca,sk[top]), sk[top]=lca;
    	sk[++top]=x;
    }
    void Build(const Vec &p)
    {
    	int n=p.size();
    	sk[top=1]=1;
    	if(p[0]==1) for(int i=1; i<n; ++i) Insert(p[i]);
    	else for(int i=0; i<n; ++i) Insert(p[i]);
    	while(--top) Add_direct(sk[top],sk[top+1]);
    }
    LL DFS(int x,int d)//虚树上有其它的点!
    {
    	LL res1=0, res2=0, s=A[x]%d?0:phi[A[x]];
    	for(int v,i=H[x]; i; i=nxt[i])
    		res1+=DFS(v=to[i],d), res2+=1ll*s*sum[v]%mod, Add(s,sum[v]);
    	sum[x]=s%mod, H[x]=0;//!
    	return (res1%mod+res2*(LL)dep[x]%mod)%mod;//子树的答案和当前点为LCA独立,别算在一块。。
    }
    void Calc(int d)
    {
    	const Vec &p=pt[d];
    //	if(!p.size()) return;
    	Build(p);
    	LL res=0, sum=0; int n=p.size();
    	for(int i=0; i<n; ++i) if(!(A[p[i]]%d)) sum+=phi[A[p[i]]];
    	sum%=mod;
    	for(int i=0; i<n; ++i) if(!(A[p[i]]%d)) res+=1ll*phi[A[p[i]]]*(sum-phi[A[p[i]]])%mod*dep[p[i]]%mod;
    	res%=mod, res-=2ll*DFS(1,d)%mod;
    	f[d]=(res%mod+mod)%mod;//+mod!
    	Enum=0;
    }
    
    int main()
    {
    	n=read();
    	for(int i=1; i<=n; ++i) A[i]=read();
    	for(int i=1; i<n; ++i) AddEdge(read(),read());
    	Pre(), preDFS1(1), preDFS2(1,1);
    	memset(H,0,sizeof H);
    	for(int i=1; i<=n; ++i) Calc(i);
    	LL ans=0;
    	for(int d=1; d<=n; ++d)
    	{
    		LL sum=0;
    		for(int i=1; i*d<=n; ++i) sum+=mu[i]*f[i*d];
    		ans+=1ll*(sum%mod)*d%mod*inv[phi[d]]%mod;
    	}
    	printf("%I64d
    ",2ll*(ans%mod)*inv[n]%mod*inv[n-1]%mod);
    
    	return 0;
    }
    
  • 相关阅读:
    【笔记】 寻址方式
    今日思考之 20200614:java 中 null 是否对 gc 有帮助?
    分布式唯一ID生成方案对比分析 笔记
    请不要再称数据库是CP或者AP——CAP的误导(短板)和它的使命
    延迟初始化中的 双重检查模式 和 延迟占位类模式 你都用对了吗?
    redis bitmap
    RabbitMQ 使用 Policies HTTP API 绑定和解绑 DLX
    spring boot 自动装配的实现原理和骚操作,不同版本实现细节,debug 到裂开......
    netty 学习笔记三:大跃进,使用 netty 实现 IM 即时通讯系统
    一道 Java 方法传值面试题——Java方法传值的值传递概念和效果 + Integer 缓存机制 + 反射修改 private final 域
  • 原文地址:https://www.cnblogs.com/SovietPower/p/9294996.html
Copyright © 2011-2022 走看看