1. 程式人生 > >[BZOJ3451] Tyvj1953 Normal 點分治+FFT

[BZOJ3451] Tyvj1953 Normal 點分治+FFT

題目連結

顯然按照點的貢獻來考慮,答案就是每個點在點分樹的期望深度之和。

深度一般可以轉化為它是幾個點的子樹。

考慮\(y\)\(x\)的點分樹的子樹中。那麼需要保證在原樹\(x\)\(y\)的路徑上的點裡面,必須是\(x\)最先被隨機到。因此這種概率為\(\frac 1{dis(x,y)}\)

所以統計一下原樹上每種長度的路徑的條數,這個東西就是點分治+FFT的套路了。

#include<bits/stdc++.h>
using namespace std;
#define fec(i,x,y) (int i=head[x],y=g[i].to;i;i=g[i].ne,y=g[i].to)
#define dbg(...) fprintf(stderr,__VA_ARGS__)
#define File(x) freopen(#x".in","r",stdin),freopen(#x".out","w",stdout)
#define isin(x,S) (((S)>>((x)-1))&1)
#define fi first
#define se second
#define pb push_back
template<typename I>inline void read(I&x){int f=0,c;while(!isdigit(c=getchar()))c=='-'?f=1:0;x=c&15;while(isdigit(c=getchar()))x=(x<<1)+(x<<3)+(c&15);f?x=-x:0;}
template<typename A,typename B>inline char SMAX(A&a,const B&b){return a<b?a=b,1:0;}
template<typename A,typename B>inline char SMIN(A&a,const B&b){return a>b?a=b,1:0;}
typedef long long ll;typedef unsigned long long ull;typedef pair<int,int>pii;

const int N=30000+7;const double PI=acos(-1);
int n,x,y,ans[N];double ans1;
struct Edge{int to,ne;}g[N<<1];int head[N],tot;
inline void Addedge(int x,int y){g[++tot].to=y;g[tot].ne=head[x];head[x]=tot;}
inline void Adde(int x,int y){Addedge(x,y),Addedge(y,x);}

struct Complex{
    double x,y;
    inline Complex(const double&x=0,const double&y=0):x(x),y(y){}
    inline Complex operator+(const Complex&a){return Complex(x+a.x,y+a.y);}
    inline Complex operator-(const Complex&a){return Complex(x-a.x,y-a.y);}
    inline Complex operator*(const Complex&a){return Complex(x*a.x-y*a.y,x*a.y+y*a.x);}
}A[N<<2],B[N<<2];

inline void FFT(Complex*a,int n,int f){
    for(int i=0,j=0;i<n;++i){
        if(i>j)swap(a[i],a[j]);
        for(int l=(n>>1);(j^=l)<l;l>>=1);
    }
    for(int i=1;i<n;i<<=1){
        Complex w(cos(PI/i),f*sin(PI/i));
        for(int j=0;j<n;j+=(i<<1)){
            Complex e(1,0);
            for(int k=0;k<i;++k,e=e*w){
                Complex x=a[j+k],y=e*a[i+j+k];
                a[j+k]=x+y,a[i+j+k]=x-y;
            }
        }
    }
    if(f>0)return;
    for(int i=0;i<n;++i)a[i].x/=n;
}
inline void Mul(int*a,int*b,int n,int m){
    int l=1;while(l<=n+m)l<<=1;
    for(int i=0;i<=n;++i)A[i]=Complex(a[i],0);
    for(int i=0;i<=m;++i)B[i]=Complex(b[i],0);
    FFT(A,l,1);FFT(B,l,1);
    for(int i=0;i<l;++i)A[i]=A[i]*B[i];
    FFT(A,l,-1);
    for(int i=0;i<=min(n+m,::n);++i)ans[i]+=(int)(A[i].x+0.5);
    for(int i=0;i<l;++i)A[i]=B[i]=Complex();
}

int rt,mima,sum,num[N],vis[N];
inline void Grt(int x,int fa=0){
    num[x]=1;int f=0;
    for fec(i,x,y)if(y!=fa&&!vis[y])Grt(y,x),num[x]+=num[y],SMAX(f,num[y]);
    SMAX(f,sum-num[x]);if(SMIN(mima,f))rt=x;
}

int f[N],s[N],fa[N],q[N],dis[N],hd,tl;
inline void BFS(int x,int f){
    q[hd=0,tl=1]=x;fa[x]=f;
    while(hd<tl){
        int x=q[++hd];dis[x]=dis[fa[x]]+1;++::f[dis[x]-1];
        for fec(i,x,y)if(y!=fa[x]&&!vis[y])q[++tl]=y,fa[y]=x;
    }
}

inline void Calc(int x){
    s[1]=1;dis[x]=1;int mxd=1;ans[s[1]]++;
    for fec(i,x,y)if(!vis[y]){
        f[0]=0;BFS(y,x);num[y]=tl;
        Mul(s,f,mxd,dis[q[tl]]-1);
        SMAX(mxd,dis[q[tl]]);
        for(int i=1;i<=dis[q[tl]];++i)s[i]+=f[i-1],f[i-1]=0;
    }
    for(int i=1;i<=mxd;++i)s[i]=0;
}

inline void Solve(int x){
    vis[x]=1;Calc(x);
    for fec(i,x,y)if(!vis[y])mima=sum=num[y],Grt(y),Solve(rt);
}

int main(){
    #ifdef hzhkk
    freopen("hkk.in","r",stdin);
    #endif
    read(n);
    for(int i=1;i<n;++i)read(x),read(y),Adde(x+1,y+1);
    sum=mima=n;Grt(1);Solve(rt);
    ans1=ans[1];
    for(int i=2;i<=n;++i)ans1+=(double)ans[i]*2/i;
    printf("%.4lf\n",ans1);
}