[BZOJ3451][Tyvj1953]Normal(点分治+FFT)

https://www.cnblogs.com/GXZlegend/p/8611948.html

#include<cmath>
#include<cstdio>
#include<algorithm>
#define rep(i,l,r) for (int i=(l); i<=(r); i++)
#define For(i,x) for (int i=h[x],k; i; i=nxt[i])
using namespace std;

const int N=100010;
const double pi=acos(-1.);
int n,u,v,S,rt,tot,cnt,mx,d[N],q[N],rev[N];
int vis[N],sz[N],f[N],h[N],to[N<<1],nxt[N<<1];
double ans,num[N];
void add(int u,int v){ to[++cnt]=v; nxt[cnt]=h[u]; h[u]=cnt; }

struct C{ double x,y; }a[N];
C operator +(const C &a,const C &b){ return (C){a.x+b.x,a.y+b.y}; }
C operator -(const C &a,const C &b){ return (C){a.x-b.x,a.y-b.y}; }
C operator *(const C &a,const C &b){ return (C){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x}; }

void find(int x,int fa){
    sz[x]=1; f[x]=0;
    For(i,x) if (!vis[k=to[i]] && k!=fa)
        find(k,x),sz[x]+=sz[k],f[x]=max(f[x],sz[k]);
    f[x]=max(f[x],S-sz[x]);
    if (f[x]<f[rt]) rt=x;
}

void get(int x,int fa){
    q[++tot]=d[x];
    For(i,x) if (!vis[k=to[i]] && k!=fa) d[k]=d[x]+1,get(k,x);
}

void DFT(C a[],int n,int f){
    for (int i=0; i<n; i++) if (i<rev[i]) swap(a[i],a[rev[i]]);
    for (int i=1; i<n; i<<=1){
        C wn=(C){cos(pi/i),f*sin(pi/i)};
        for (int p=i<<1,j=0; j<n; j+=p){
            C w=(C){1,0};
            for (int k=0; k<i; k++,w=w*wn){
                C x=a[j+k],y=w*a[i+j+k]; a[j+k]=x+y; a[i+j+k]=x-y;
            }
        }
    }
    if (f==-1) for (int i=0; i<n; i++) a[i].x/=n;
}

void calc(int fl){
    int n=1,L=0; mx=0;
    rep(i,1,tot) mx=max(mx,q[i]);
    for (; n<=mx<<1; n<<=1) L++;
    for (int i=0; i<n; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
    for (int i=0; i<n; i++) a[i].x=a[i].y=0;
    rep(i,1,tot) a[q[i]].x++;
    DFT(a,n,1);
    for (int i=0; i<n; i++) a[i]=a[i]*a[i];
    DFT(a,n,-1);
    rep(i,0,2*mx) num[i]+=fl*(int)(a[i].x+0.1);
}

void work(int x){
    vis[x]=1; d[x]=tot=0; get(x,0); calc(1);
    For(i,x) if (!vis[k=to[i]]){
        d[k]=1; tot=0; get(k,0); calc(-1);
        S=sz[k]; rt=0; find(k,0); work(rt);
    }
}

int main(){
    freopen("bzoj3451.in","r",stdin);
    freopen("bzoj3451.out","w",stdout);
    scanf("%d",&n);
    rep(i,2,n) scanf("%d%d",&u,&v),add(u+1,v+1),add(v+1,u+1);
    S=f[0]=n; find(1,0); work(rt);
    rep(i,1,n) ans+=num[i-1]/i;
    printf("%.4lf
",ans);
    return 0;
}