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

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

 1 #include<cmath>
 2 #include<cstdio>
 3 #include<algorithm>
 4 #define rep(i,l,r) for (int i=(l); i<=(r); i++)
 5 #define For(i,x) for (int i=h[x],k; i; i=nxt[i])
 6 using namespace std;
 7 
 8 const int N=100010;
 9 const double pi=acos(-1.);
10 int n,u,v,S,rt,tot,cnt,mx,d[N],q[N],rev[N];
11 int vis[N],sz[N],f[N],h[N],to[N<<1],nxt[N<<1];
12 double ans,num[N];
13 void add(int u,int v){ to[++cnt]=v; nxt[cnt]=h[u]; h[u]=cnt; }
14 
15 struct C{ double x,y; }a[N];
16 C operator +(const C &a,const C &b){ return (C){a.x+b.x,a.y+b.y}; }
17 C operator -(const C &a,const C &b){ return (C){a.x-b.x,a.y-b.y}; }
18 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}; }
19 
20 void find(int x,int fa){
21     sz[x]=1; f[x]=0;
22     For(i,x) if (!vis[k=to[i]] && k!=fa)
23         find(k,x),sz[x]+=sz[k],f[x]=max(f[x],sz[k]);
24     f[x]=max(f[x],S-sz[x]);
25     if (f[x]<f[rt]) rt=x;
26 }
27 
28 void get(int x,int fa){
29     q[++tot]=d[x];
30     For(i,x) if (!vis[k=to[i]] && k!=fa) d[k]=d[x]+1,get(k,x);
31 }
32 
33 void DFT(C a[],int n,int f){
34     for (int i=0; i<n; i++) if (i<rev[i]) swap(a[i],a[rev[i]]);
35     for (int i=1; i<n; i<<=1){
36         C wn=(C){cos(pi/i),f*sin(pi/i)};
37         for (int p=i<<1,j=0; j<n; j+=p){
38             C w=(C){1,0};
39             for (int k=0; k<i; k++,w=w*wn){
40                 C x=a[j+k],y=w*a[i+j+k]; a[j+k]=x+y; a[i+j+k]=x-y;
41             }
42         }
43     }
44     if (f==-1) for (int i=0; i<n; i++) a[i].x/=n;
45 }
46 
47 void calc(int fl){
48     int n=1,L=0; mx=0;
49     rep(i,1,tot) mx=max(mx,q[i]);
50     for (; n<=mx<<1; n<<=1) L++;
51     for (int i=0; i<n; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
52     for (int i=0; i<n; i++) a[i].x=a[i].y=0;
53     rep(i,1,tot) a[q[i]].x++;
54     DFT(a,n,1);
55     for (int i=0; i<n; i++) a[i]=a[i]*a[i];
56     DFT(a,n,-1);
57     rep(i,0,2*mx) num[i]+=fl*(int)(a[i].x+0.1);
58 }
59 
60 void work(int x){
61     vis[x]=1; d[x]=tot=0; get(x,0); calc(1);
62     For(i,x) if (!vis[k=to[i]]){
63         d[k]=1; tot=0; get(k,0); calc(-1);
64         S=sz[k]; rt=0; find(k,0); work(rt);
65     }
66 }
67 
68 int main(){
69     freopen("bzoj3451.in","r",stdin);
70     freopen("bzoj3451.out","w",stdout);
71     scanf("%d",&n);
72     rep(i,2,n) scanf("%d%d",&u,&v),add(u+1,v+1),add(v+1,u+1);
73     S=f[0]=n; find(1,0); work(rt);
74     rep(i,1,n) ans+=num[i-1]/i;
75     printf("%.4lf\n",ans);
76     return 0;
77 }

 

转载于:https://www.cnblogs.com/HocRiser/p/10513358.html

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值