bzoj3451 Normal （点分治+FFT）

根据期望的线性性，考虑每个点的点分树子树期望大小

j 会给 i 贡献 1 的大小当且仅当 j 在 i 的子树内

所以就是这个式子 ∑i∑jp[i][j]  p[i][j]表示 j 在 i 子树内的概率

j 在 i 的点分树子树内，相当于（i,j）这个路径上，i 是第一个被选择的分治重心

而这个概率就是 p[i][j] = 1/(dis(i,j)+1)

统计出所有距离为 i 的点对数量，这个可以用点分治来求，统计的时候相当于在做卷积，可以用 fft 优化

时间复杂度O(nlog^2n)

#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<queue>
#include<stack>
#include<cmath>
using namespace std;
typedef long long ll;

const int maxn = 30010;
const int INF = 1e9+7;
const double Pi = acos(-1.0);

int n,m;
int maxson,sum,rt;
int ans[maxn];
double res=0;

int h[maxn],size;
struct E{
    int to,next;
}e[maxn<<1];
void add(int u,int v){
    e[++size].to=v;
    e[size].next=h[u];
    h[u]=size;
}

int rev[maxn],lim=1,l=0;

struct Complex{
    double x,y;
    Complex(double xx=0,double yy=0){ x=xx,y=yy; }
}A[maxn],B[maxn];
Complex operator + (Complex a,Complex b){ return Complex(a.x+b.x,a.y+b.y); }
Complex operator - (Complex a,Complex b){ return Complex(a.x-b.x,a.y-b.y); }
Complex operator * (Complex a,Complex b){ return Complex(a.x*b.x-a.y*b.y,a.y*b.x+a.x*b.y); }

void fft(Complex *A,int type){
    for(int i=0;i<lim;i++) if(i<rev[i]) swap(A[i],A[rev[i]]);
    for(int i=1;i<lim;i<<=1){
        Complex Wn(cos(Pi/i),type*sin(Pi/i));
        for(int j=0,p=i<<1;j<lim;j+=p){
            Complex w(1,0);
            for(int k=0;k<i;k++,w=w*Wn){
                Complex x=A[j+k],y=w*A[i+j+k];
                A[j+k]=x+y;
                A[i+j+k]=x-y;
            }
        }
    }
//    if(type==-1) for(int i=0;i<lim;i++) A[i].x/=lim;
}

int vis[maxn];
int sz[maxn],son[maxn];
void getroot(int u,int par){
    sz[u]=1; son[u]=0;
    for(int i=h[u];i!=-1;i=e[i].next){
        int v=e[i].to;
        if(v==par||vis[v]) continue;
        getroot(v,u);
        sz[u]+=sz[v];
        son[u]=max(son[u],sz[v]);
    }
    son[u]=max(son[u],sum-sz[u]);
    if(maxson>son[u]){
        maxson=son[u];
        rt=u;
    }
}

int t[maxn],s[maxn],tmp[maxn],mxd,md;

void dfs(int u,int par,int dep){
    t[dep]++; md=max(md,dep);
    for(int i=h[u];i!=-1;i=e[i].next){
        int v=e[i].to;
        if(vis[v]||v==par) continue;
        dfs(v,u,dep+1);
    }
}

void mul(int *a,int *b,int _n,int _m,int *c){
    lim=1,l=0;
    while(lim<=_n+_m) lim<<=1,l++;
    for(int i=0;i<lim;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
    for(int i=0;i<=_n;i++) A[i].x=1.0*a[i];
    for(int i=0;i<=_m;i++) B[i].x=1.0*b[i];
    
    fft(A,1); fft(B,1);
    for(int i=0;i<lim;i++) A[i]=A[i]*B[i];
    fft(A,-1);
    
    for(int i=0;i<=_n+_m;i++){ c[i]=(int)(A[i].x/lim+0.5); }
    for(int i=0;i<lim;i++) A[i]=B[i]=(Complex){0,0};
}

void divide(int u){
    vis[u]=1; int mxd=0;
    s[0]=1;
    for(int i=h[u];i!=-1;i=e[i].next){
        int v=e[i].to;
        if(vis[v]) continue;
        md=0; dfs(v,u,1);
        
        mul(s,t,mxd,md,tmp);
//        for(int j=0;j<=md;j++) printf("%d ",tmp[j]); printf("\n");
        for(int j=0;j<=mxd+md;j++) ans[j]+=tmp[j],tmp[j]=0;
        for(int j=0;j<=md;j++) s[j]+=t[j],t[j]=0; // 合并已经遍历过的子树中的信息
        mxd=max(mxd,md);
    }
//    printf("zhixing\n");
    for(int i=0;i<=mxd;i++) s[i]=0;
    for(int i=h[u];i!=-1;i=e[i].next){
        int v=e[i].to;
        if(vis[v]) continue;
        rt=0,sum=sz[v],maxson=INF;
        getroot(v,u);
        divide(rt);
    }
}

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

int main(){
    memset(h,-1,sizeof(h));
    n=read();
    int u,v;
    for(int i=1;i<n;i++){
        u=read(),v=read();
        add(u,v),add(v,u);
    }
    
    sum=n;
    maxson=INF;
    rt=0;
    getroot(1,0);

    divide(rt);
    
//    for(int i=1;i<=n;i++) printf("%d ",ans[i]); printf("\n"); 
    
    for(int i=1;i<=n;i++) res+=1.0*ans[i]/(i+1);
    res*=2; // 路径是相互的
    res+=n; // 自己对自己的贡献 
    printf("%.4lf\n",res);

    return 0;
}