Surprise Me![莫比乌斯反演+欧拉函数+虚树]
最近在学反演,想起来考过一道难得不纯推式子的莫反,特地抠下来水一篇博。
首先忽略方案数,最后统一除去。
对欧拉函数 \(\varphi(x)\) ,有 $$\varphi(x\times y)=\frac{\gcd(x,y)\times\varphi(x)\times\varphi(y)}{\varphi(\gcd(x,y))}$$
可以从 \(\varphi(x)=x\times\prod_{p|x}\frac{p-1}{p}\) 的角度理解。两个欧拉函数相乘后多计算的部分即 \(\gcd\) 的质因子。
于是推式子。令 \(a_{loc_i}=i\) ,即编号为 \(loc_i\) 的节点数值为 \(i\) ,有
莫比乌斯反演,有 \([\gcd(i,j)=1]=\sum_{d|\gcd(i,j)}\mu(d)\) ,于是有
设 \(t=gd\) ,那么
到这里就可以了。 \(\sum_{g|t}\frac{g\mu(\frac{t}{g})}{\varphi(g)}\) 是关于 \(t\) 的函数,可以线性筛出 \(\mu\) 与 \(\varphi\) 后枚举倍数 \(O(n\ln n)\) 得出。考虑后面的东西怎么求。
与上面相似的思想,同样枚举 \(t\) ,每次考虑 \(t\) 的倍数的数值对答案的贡献。这样复杂度同样是 \(O(n\ln n)\) 。
计算贡献时,因为每次考虑的节点不多,可以建出虚树后 \(DP\) 。令每个点的权值 \(w_i=\varphi(a_i)[t|i]\) ,即关键点的点权为它的数值的欧拉函数,其他点权为 \(0\) ,那么我们要求的就是 \(\sum_i\sum_jw_iw_j\operatorname{dist}(i,j)\) 。
继续拆式子,设 \(v\) 是 \(s\) 的儿子,现将 \(v\) 合并至 \(s\) ,设 \(g_s\) 为 \(s\) 子树中的答案,那么 \(g_s=g_s+g_v+delta\) ,
这样就把 \(x,s\) 和 \(y,v\) 放到了一起。设
- \(f_s=\sum_{x\in s}w_x\operatorname{dist}(x,s)\)
- \(d_s=\sum_{x\in s}w_x\)
进行辅助,就能完成转移。
\(code:\)
T1
#include<bits/stdc++.h>
using namespace std;
namespace IO{
typedef long long LL;
#define int LL
int read(){
int x=0,f=0; char ch=getchar();
while(ch<'0'||ch>'9'){ f|=(ch=='-'); ch=getchar(); }
while(ch>='0'&&ch<='9'){ x=(x<<1)+(x<<3)+(ch^48); ch=getchar(); }
return f?-x:x;
} char output[50];
void write(LL x,char sp){
int len=0;
if(x<0) x=-x, putchar('-');
do{ output[len++]=x%10+'0'; x/=10; }while(x);
for(int i=len-1;~i;i--) putchar(output[i]); putchar(sp);
}
void ckmin(int &x,int y){ x=x<y?x:y; }
void ckmax(int &x,int y){ x=x>y?x:y; }
} using namespace IO;
const int NN=200010,mod=1e9+7;
int n,ans,p[NN],loc[NN];
vector<int>e[NN];
void add(int a,int b){
e[a].push_back(b);
e[b].push_back(a);
}
namespace Mathematics{
int cnt,mu[NN],pri[NN],phi[NN],bas[NN];
bool vis[NN];
int qpow(int a,int b,int res=1){
for(;b;b>>=1,a=a*a%mod)
if(b&1) res=res*a%mod;
return res;
}
void get_functions(){
mu[1]=phi[1]=1;
for(int i=2;i<=n;i++){
if(!vis[i]) mu[i]=-1, phi[i]=i-1, pri[++cnt]=i;
for(int j=1;j<=cnt&&pri[j]*i<=n;j++){
int now=pri[j]*i;
vis[now]=1;
if(i%pri[j]==0){
mu[now]=0; phi[now]=phi[i]*pri[j];
break;
}
mu[now]=-mu[i]; phi[now]=phi[i]*phi[pri[j]];
}
}
for(int i=1;i<=n;i++) for(int j=i;j<=n;j+=i)
(bas[j]+=mod+mu[j/i]*i*qpow(phi[i],mod-2)%mod)%=mod;
}
} using namespace Mathematics;
namespace RMQforLCA{
int idx,fa[NN],lg[NN<<1],dep[NN],dfn[NN],seq[NN<<1][21];
int LCA(int x,int y){
x=dfn[x]; y=dfn[y];
if(x>y) swap(x,y);
int t=lg[y-x+1],u=seq[x][t],v=seq[y-(1<<t)+1][t];
return dep[u]<dep[v]?u:v;
}
void RMQ_dfs(int s,int f){
dfn[s]=++idx; fa[s]=f; seq[idx][0]=s; dep[s]=dep[f]+1;
for(int v:e[s]) if(v!=f)
RMQ_dfs(v,s), seq[++idx][0]=s;
}
void build_RMQ(){
for(int i=2;i<=idx;i++) lg[i]=lg[i>>1]+1;
for(int i=1;i<=20;i++)
for(int u,v,j=1;j<=idx-(1<<i)+1;j++){
u=seq[j][i-1]; v=seq[j+(1<<i-1)][i-1];
seq[j][i]=dep[u]<dep[v]?u:v;
}
}
} using namespace RMQforLCA;
namespace Virtual_Tree{
int w[NN],s[NN],f[NN],g[NN];
int top,stk[NN];
vector<int>ve[NN],key;
bool cmp(int a,int b){ return dfn[a]<dfn[b]; }
void vadd(int a,int b){ ve[a].push_back(b); }
void insert(int x){
if(top==1) return stk[++top]=x,void();
int lca=LCA(x,stk[top]);
if(lca==stk[top]) return stk[++top]=x,void();
while(top>1&&dfn[stk[top-1]]>=dfn[lca]) vadd(stk[top-1],stk[top]), --top;
if(lca!=stk[top]) vadd(lca,stk[top]), stk[top]=lca;
stk[++top]=x;
}
void build(int t){
key.clear(); stk[top=1]=1;
for(int i=t;i<=n;i+=t){
w[loc[i]]=phi[i];
if(loc[i]^1) key.push_back(loc[i]);
}
sort(key.begin(),key.end(),cmp);
for(int x:key) insert(x);
while(top>1) vadd(stk[top-1],stk[top]), --top;
}
void dfs(int u){
s[u]=w[u]; f[u]=g[u]=0;
for(int v:ve[u]) if(v!=fa[u]){
int len=dep[v]-dep[u];
dfs(v);
(g[u]+=g[v]+f[u]*s[v]+s[u]*s[v]%mod*len+f[v]*s[u])%=mod;
(f[u]+=s[v]*len+f[v])%=mod;
(s[u]+=s[v])%=mod;
}
ve[u].clear(); w[u]=0;
}
} using namespace Virtual_Tree;
signed main(){
n=read(); get_functions();
for(int i=1;i<=n;i++) p[i]=read(), loc[p[i]]=i;
for(int i=1;i<n;i++) add(read(),read());
RMQ_dfs(1,0); build_RMQ();
for(int i=1;i<=n;i++){
build(i); dfs(1);
(ans+=bas[i]*g[1]*2%mod)%=mod;
}
ans=ans*qpow(n,mod-2)%mod*qpow(n-1,mod-2)%mod;
write(ans,'\n');
return 0;
}

浙公网安备 33010602011771号