bzoj3451/Tyvj1953：Normal（点分治+FFT）

题面
题意：给你一棵树，XJB选点来做点分治，问期望的复杂度。

根据我对期望的粗鄙理解，期望就是个积分，所以它满足积分加减等线性运算。故这里可以考虑对每个点求期望复杂度。

对于点x，它产生的期望就是x在点分树中的期望深度。

再次运用期望的线性性质，我们可以对于每个y，求出x在y子树中的概率，由此推出期望。

经过小小的分析，发现x在y的子树中，仅当在x-y的路径上，y被第一个选中来分治，则概率为$\frac{1}{dis(x,y)}$，也就是y能为x带来的期望。

然后原问题的答案为

$$\sum_{i=1}^n\sum_{j=1}^{n}\frac{1}{dis(i,j)}$$

若dis(i,j)不在分母里，就是一个果的点分治了。

我们换下枚举的顺序，设f[x]为距离为x的点对数，变成$\sum_{x=1}^{n}f[x]$，就可以算答案了。

同样考虑点分治，找到重心x，考虑过x的路径。
对于每个连通块。求出它们每个点到x的距离，设g[i]为到x距离为i的点数。我们发现两个g能对答案的贡献是个卷积的形式，用fft优化即可。

期望复杂度O($nlog^2n$)，但理论上说，不加优化的话，是会被菊花图卡的。
比如说这样

用一条长链把fft的次数界撑大了，然后用大大的次数界不断地做fft。
应该先做深度小的连通块，再做深度大的。

#include <iostream>
#include <fstream>
#include <algorithm>
#include <cmath>
#include <ctime>
#include <cstdio>
#include <cstdlib>
#include <cstring>

using namespace std;
#define mmst(a, b) memset(a, b, sizeof(a))
#define mmcp(a, b) memcpy(a, b, sizeof(b))

typedef long long LL;

const int N=60060,oo=1e9+7;
const LL g=3,p=1004535809;

int n,rev[N];
int nn;
int to[N],nex[N],head[N],cnt;
int root,sum,siz[N],d[N],dep[N],len,slen;
LL a[N],b[N],c[N],ans[N];
double res=0.0;
bool vis[N];

void read(int &hy)
{
    hy=0;
    char cc=getchar();
    while(cc<'0'||cc>'9')
    cc=getchar();
    while(cc>='0'&&cc<='9')
    {
        hy=(hy<<3)+(hy<<1)+cc-'0';
        cc=getchar();
    }
}

LL cheng(LL a,LL b)
{
    LL res=1ll;
    for(;b;b>>=1,a=a*a%p)
    if(b&1)
    res=res*a%p;
    return res;
}

void init(int lim)
{
    int k=-1;
    n=1;
    while(n<lim)
    n<<=1,k++;
    for(int i=0;i<n;i++)
    rev[i]=(rev[i>>1] >> 1) | ((i&1)<<k);
}

void ntt(LL *a,int ops)
{
    for(int i=0;i<n;i++)
    if(i<rev[i])
    swap(a[i],a[rev[i]]);
    for(int l=2;l<=n;l<<=1)
    {
        int m=l>>1;
        LL wn;
        if(ops)
        wn=cheng(g,(p-1)/l);
        else
        wn=cheng(g,p-1-(p-1)/l);
        for(int i=0;i<n;i+=l)
        {
            LL w=1ll;
            for(int k=0;k<m;k++)
            {
                LL t=a[i+k+m]*w%p;
                a[i+k+m]=(a[i+k]-t+p)%p;
                a[i+k]=(a[i+k]+t)%p;
                w=w*wn%p;
            }
        }
    }
    if(!ops)
    {
        LL Inv=cheng(n,p-2);
        for(int i=0;i<n;i++)
        a[i]=a[i]*Inv%p;
    }
}

void add(int u,int v)
{
    to[++cnt]=v;
    nex[cnt]=head[u];
    head[u]=cnt;
}

void dfsRoot(int x,int fa)
{
    d[x]=0;
    siz[x]=1;
    for(int h=head[x];h;h=nex[h])
    if(!vis[to[h]]&&to[h]!=fa)
    {
        dfsRoot(to[h],x);
        siz[x]+=siz[to[h]];
        d[x]=max(d[x],siz[to[h]]);
    }
    d[x]=max(d[x],sum-siz[x]);
    if(d[x]<d[root])
    root=x;
}

void dfsLen(int x,int fa)
{
    len=max(len,dep[x]);
    a[dep[x]]++;
    for(int h=head[x];h;h=nex[h])
    if(!vis[to[h]]&&to[h]!=fa)
    {
        dep[to[h]]=dep[x]+1;    
        dfsLen(to[h],x);
    }
}

void dfsSol(int x)
{
    vis[x]=1;
    slen=1;
    b[0]=1;
    for(int h=head[x];h;h=nex[h])
    if(!vis[to[h]])
    {
        len=0;
        dep[to[h]]=1;
        dfsLen(to[h],0);

        init(len+slen+1);
        for(int i=0;i<n;i++)
        c[i]=b[i];

        for(int i=0;i<=len;i++)
        b[i]+=a[i];

        ntt(c,1);
        ntt(a,1);

        for(int i=0;i<n;i++)
        c[i]=c[i]*a[i]%p;

        ntt(c,0);
        for(int i=0;i<n;i++)
        ans[i]+=c[i];

        for(int i=0;i<n;i++)
        a[i]=0;
        slen=max(slen,len);
    }
    for(int i=0;i<=slen;i++)
    b[i]=0;
    for(int h=head[x];h;h=nex[h])
    if(!vis[to[h]])
    {
        root=0;
        sum=siz[to[h]];
        dfsRoot(to[h],0);
        dfsSol(root);
    }
}

int main()
{
    cin>>nn;

    for(int i=1;i<nn;i++)
    {
        int u,v;
        read(u);
        read(v);
        u++;
        v++;
        add(u,v);
        add(v,u);
    }
    sum=nn;
    root=0;
    d[0]=oo;
    dfsRoot(1,0);
    dfsSol(root);

    for(int i=1;i<nn;i++)
    res+=2.0*ans[i]/(i+1);
    res+=nn;
    printf("%.4lf\n",res);
}