Codeforces 809E：Surprise me! （莫比乌斯反演+虚树）

题目传送门：http://codeforces.com/contest/809/problem/E

---

题目分析：一道套路到极点的题目。

公式推导直接见https://blog.sengxian.com/solutions/cf-809e，我懒得打了QAQ。

最后推出：

$$ans=\sum_{T=1}^{n}\sum_{d|T}\frac{d\mu(\frac{T}{d})}{\phi(d)}\sum_{d|a_i}\sum_{d|a_j}\phi(a_i)\phi(a_j)dis(i,j)$$

把第二个$\sum$后面的玩意儿记为$G(T)$，这个可以用$O(n\ln(n))$预处理出来。然后对于所有满足$d|a_i$的点i建虚树，枚举虚树中的任意一条边，用这条边的长度乘以它左右两边的$\phi$值之积贡献答案。由于a是1~n的一个排列，虚树的总大小是$n\ln(n)$的。用ST表预处理可以做到$O(1)$求LCA。如果对所有满足条件的点i按时间戳排序，会使时间复杂度多个$\log$。可以一开始先对所有点排序，然后按顺序插入vector中。由于要枚举因数，时间是$\sum_{i=1}^{n}\sqrt{i}$的，大概是$6*10^7$（其实和排序的时间差不了多少不过我还是这样写了）。

一开始以为虚树的边边权是1(naive!)，然后WA了一次。后来RE了，我还以为是vector的问题（有的时候STL的容器太大会蜜汁RE），后来发现是自己空间池开小了，我没有考虑建虚树时的加边。

所以我做这题就是为了练习写虚树的，这大概是最后一道虚树的题了

然后我把CODE中虚树的模板抽了出来：

tail=1;
sak[1]=tree[1];
for (int i=2; i<=Node[d].size(); i++)
{
    int x=tree[i],last=0,p=Lca(sak[tail],x);
    while ( dep[ sak[tail] ]>dep[p] && tail ) last=sak[tail--];
    if (sak[tail]!=p) Fa[p]=sak[tail],sak[++tail]=tree[++cnt]=p;
    if (last) Fa[last]=p;
    Fa[x]=p;
    sak[++tail]=x; //!!!
}
int root=sak[1];
Fa[root]=0;

由于这代码是自己根据对虚树的理解YY出来的，所以十分简(chou)洁(lou)，纯属给自己有空的时候看看，不喜勿喷。

这题本来是计划1h过的，结果1h的时候刚刚敲完代码，又过了25min才调过，代码能力还是要加强。

---

CODE：

#include<iostream>
#include<string>
#include<cstring>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<stdio.h>
#include<algorithm>
#include<vector>
using namespace std;

const int maxn=200100;
const int maxl=20;
const long long M=1000000007;
typedef long long LL;

struct edge
{
    int obj;
    edge *Next;
} e[maxn<<1];
edge *head[maxn];
int cur=0;

vector <int> Node[maxn];

int Rank[maxn];
int dep[maxn];
int dfn[maxn];
int Time=0;

int Lg[maxn<<1];
int dfsx[maxn<<1];
int ST[maxn<<1][maxl];

int tree[maxn];
int cnt;
int sak[maxn];
int tail;

LL Size[maxn];
int Fa[maxn];

int phi[maxn];
int miu[maxn];

LL rev[maxn];
LL G[maxn];

bool vis[maxn];
int prime[maxn];
int num=0;

int a[maxn];
int n;
LL ans=0;

void Add(int x,int y)
{
    cur++;
    e[cur].obj=y;
    e[cur].Next=head[x];
    head[x]=e+cur;
}

void Linear_shaker()
{
    phi[1]=miu[1]=1;
    for (int i=2; i<=n; i++)
    {
        if (!vis[i]) phi[i]=i-1,miu[i]=-1,prime[++num]=i;
        for (int j=1; j<=num && i*prime[j]<=n; j++)
        {
            int k=i*prime[j];
            vis[k]=true;
            if (i%prime[j])
            {
                phi[k]=phi[i]*(prime[j]-1);
                miu[k]=-miu[i];
            }
            else
            {
                phi[k]=phi[i]*prime[j];
                miu[k]=0;
                break;
            }
        }
    }

    rev[1]=1;
    for (LL i=2; i<=n; i++)
    {
        LL x=M/i,y=M%i;
        rev[i]=x*rev[y]%M;
        rev[i]=M-rev[i];
    }

    for (int i=1; i<=n; i++) if (miu[i]==-1) miu[i]=M-1;
    for (int i=1; i<=n; i++)
        for (int j=1; i*j<=n; j++)
            G[i*j]=(G[i*j]+ (long long)i*miu[j]%M*rev[ phi[i] ]%M )%M;
}

void Dfs(int node,int fa)
{
    dfn[node]=++Time;
    dfsx[Time]=node;
    for (edge *p=head[node]; p; p=p->Next)
    {
        int son=p->obj;
        if (son==fa) continue;
        dep[son]=dep[node]+1;
        Dfs(son,node);
        dfsx[++Time]=node;
    }
}

bool Comp(int x,int y)
{
    return dfn[x]<dfn[y];
}

int Lca(int x,int y)
{
    x=dfn[x],y=dfn[y];
    if (x>y) swap(x,y);
    int lg=Lg[y-x];
    y=y-(1<<lg)+1;
    x=ST[x][lg];
    y=ST[y][lg];
    if (dep[x]<dep[y]) return x;
    return y;
}

void Calc(int node)
{
    for (edge *p=head[node]; p; p=p->Next)
    {
        int son=p->obj;
        Calc(son);
        Size[node]=(Size[node]+Size[son])%M;
    }
}

int main()
{
    freopen("c.in","r",stdin);
    freopen("c.out","w",stdout);

    scanf("%d",&n);
    for (int i=1; i<=n; i++) scanf("%d",&a[i]),head[i]=NULL;
    for (int i=1; i<n; i++)
    {
        int x,y;
        scanf("%d%d",&x,&y);
        Add(x,y);
        Add(y,x);
    }

    Linear_shaker();

    dep[1]=1;
    Dfs(1,1);

    for (int i=1; i<=Time; i++) ST[i][0]=dfsx[i];
    for (int j=1; j<maxl; j++)
        for (int i=1; i<=Time; i++)
        {
            int mid=min(i+(1<<(j-1)),Time);
            if (dep[ ST[i][j-1] ]<dep[ ST[mid][j-1] ]) ST[i][j]=ST[i][j-1];
            else ST[i][j]=ST[mid][j-1];
        }
    Lg[1]=0;
    for (int i=2; i<=Time; i++) Lg[i]=Lg[i>>1]+1;

    for (int i=1; i<=n; i++) Rank[i]=i;
    sort(Rank+1,Rank+n+1,Comp);
    for (int i=1; i<=n; i++)
    {
        int node=Rank[i];
        int max_d=(int)floor( sqrt( (double)a[node] )+0.001 );
        for (int d=1; d<=max_d; d++)
            if (a[node]%d==0)
            {
                Node[d].push_back(node);
                if (d!=a[node]/d) Node[ a[node]/d ].push_back(node);
            }
    }

    for (int d=1; d<=n; d++)
    {
        cnt=0;
        for (int i=0; i<Node[d].size(); i++) tree[++cnt]=Node[d][i];

        tail=1;
        sak[1]=tree[1];
        for (int i=2; i<=Node[d].size(); i++)
        {
            int x=tree[i],last=0,p=Lca(sak[tail],x);
            while ( dep[ sak[tail] ]>dep[p] && tail ) last=sak[tail--];
            if (sak[tail]!=p) Fa[p]=sak[tail],sak[++tail]=tree[++cnt]=p;
            if (last) Fa[last]=p;
            Fa[x]=p;
            sak[++tail]=x; //!!!
        }
        int root=sak[1];
        Fa[root]=0;

        cur=0;
        for (int i=1; i<=cnt; i++) head[ tree[i] ]=NULL;
        for (int i=1; i<=cnt; i++)
        {
            int x=tree[i];
            if (Fa[x]) Add(Fa[x],x);
            if (i<=Node[d].size()) Size[x]=phi[ a[x] ];
            else Size[x]=0;
        }

        Calc(root);
        LL tot=Size[root],temp=0;
        for (int i=1; i<=cnt; i++)
        {
            int x=tree[i];
            if (x==root) continue;
            temp=(temp+ Size[x]*(tot-Size[x]+M)%M*(dep[x]-dep[ Fa[x] ])%M )%M; //!!!
        }
        temp=temp*G[d]%M;
        ans=(ans+temp)%M;
    }

    ans=ans*rev[n]%M;
    ans=ans*rev[n-1]%M;
    ans=ans*2LL%M;
    printf("%I64d\n",ans);

    return 0;
}