bzoj 3451: Tyvj1953 Normal 点分治+fft

题意

某天WJMZBMR学习了一个神奇的算法：树的点分治！
这个算法的核心是这样的：
消耗时间=0
Solve(树 a)
消耗时间 += a 的 大小
如果 a 中 只有 1 个点
退出
否则在a中选一个点x，在a中删除点x
那么a变成了几个小一点的树，对每个小树递归调用Solve
我们注意到的这个算法的时间复杂度跟选择的点x是密切相关的。
如果x是树的重心，那么时间复杂度就是O(nlogn)
但是由于WJMZBMR比较傻逼，他决定随机在a中选择一个点作为x！
Sevenkplus告诉他这样做的最坏复杂度是O(n^2)
但是WJMZBMR就是不信>_<。。。
于是Sevenkplus花了几分钟写了一个程序证明了这一点。。。你也试试看吧^_^
现在给你一颗树，你能告诉WJMZBMR他的傻逼算法需要的期望消耗时间吗？（消耗时间按在Solve里面的那个为标准）
n<=30000

分析

考虑点x对点y的贡献，也就是y是x点分树祖先的期望是多少。分析一下不难发现只要x到y路径上y是第一个被选到的，那么y就是x点分树的祖先，所以期望就是1/dis(x,y)，其中dis(x,y)表示x到y路径上的点数。
那么我们要求的就是$\sum_{i=1}^n\sum_{j=1}^n\frac{1}{dis(i,j)}$。
考虑用点分治来做，然后统计某种长度的路径数量可以用fft。
这里点分治有两种写法，一种是每遍历完一棵子树就把该子树的贡献合并上来。这样的话每次都要做一遍fft，复杂度大概是O(n^2log^2)。
还有一种写法是先搜分治中心所在的整棵树，算完贡献后再搜分治中心的每棵子树，然后减掉子树内的贡献。
这样复杂度就可以到达O(nlog^2)了。

代码

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;

typedef long double db;

const int N=30005;
const double pi=acos(-1);

int n,cnt,last[N],size[N],f[N],root,tim[N],tim_now,mx_dep,dis[N*4],ans[N*2],rev[N*4],sum;
bool vis[N];
struct edge{int to,next;}e[N*2];
struct com
{
    double x,y;
    com operator + (const com &a) {return (com){x+a.x,y+a.y};}
    com operator - (const com &a) {return (com){x-a.x,y-a.y};}
    com operator * (const com &a) {return (com){x*a.x-y*a.y,a.x*y+a.y*x};}
    com operator / (const double &a) {return (com){x/a,y/a};}
};
com a[N*4];

int read()
{
    int x=0,f=1;char ch=getchar();
    while (ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while (ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
    return x*f;
}

void addedge(int u,int v)
{
    e[++cnt].to=v;e[cnt].next=last[u];last[u]=cnt;
    e[++cnt].to=u;e[cnt].next=last[v];last[v]=cnt;
}

void get_root(int x,int fa)
{
    size[x]=1;f[x]=0;
    for (int i=last[x];i;i=e[i].next)
    {
        if (e[i].to==fa||vis[e[i].to]) continue;
        get_root(e[i].to,x);
        size[x]+=size[e[i].to];
        f[x]=max(f[x],size[e[i].to]);
    }
    f[x]=max(f[x],sum-size[x]);
    if (f[x]<f[root]||!root) root=x;
}

void dfs(int x,int fa,int d)
{
    if (tim[d]!=tim_now) tim[d]=tim_now,dis[d]=0;
    dis[d]++;mx_dep=max(mx_dep,d);
    size[x]=1;
    for (int i=last[x];i;i=e[i].next)
    {
        if (e[i].to==fa||vis[e[i].to]) continue;
        dfs(e[i].to,x,d+1);
        size[x]+=size[e[i].to];
    }
}

void fft(com *a,int f,int L)
{
    for (int i=0;i<L;i++) if (i<rev[i]) swap(a[i],a[rev[i]]);
    for (int i=1;i<L;i<<=1)
    {
        com wn=(com){cos(pi/i),f*sin(pi/i)};
        for (int j=0;j<L;j+=(i<<1))
        {
            com w=(com){1,0};
            for (int k=0;k<i;k++)
            {
                com u=a[j+k],v=w*a[j+k+i];
                a[j+k]=u+v;a[j+k+i]=u-v;
                w=w*wn;
            }
        }
    }
    if (f==-1) for (int i=0;i<L;i++) a[i]=a[i]/L;
}

void calc(int op)
{
    int L=1,lg=0;
    for (;L<=mx_dep*2;L<<=1,lg++);
    for (int i=0;i<L;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(lg-1));
    for (int i=mx_dep+1;i<L;i++) dis[i]=0;
    for (int i=0;i<L;i++) a[i]=(com){dis[i],0.0};
    fft(a,1,L);
    for (int i=0;i<L;i++) a[i]=a[i]*a[i];
    fft(a,-1,L);
    for (int i=1;i<=mx_dep*2;i++) ans[i]+=op*(int)(a[i].x+0.1);
}

void solve(int x)
{
    mx_dep=0;tim_now++;dfs(x,0,0);calc(1);
    vis[x]=1;
    for (int i=last[x];i;i=e[i].next)
    {
        if (vis[e[i].to]) continue;
        mx_dep=0;tim_now++;dis[0]=0;
        dfs(e[i].to,x,1);calc(-1);
        root=0;sum=size[e[i].to];get_root(e[i].to,x);
        solve(root);
    }
}

int main()
{
    n=read();
    for (int i=1;i<n;i++)
    {
        int x=read()+1,y=read()+1;
        addedge(x,y);
    }
    root=0;sum=n;get_root(1,0);
    solve(root);
    db anss=n;
    for (int i=1;i<=n*2;i++) if (ans[i]) anss+=(db)1/(i+1)*ans[i];
    printf("%.4lf",(double)anss);
    return 0;
}