zoukankan      html  css  js  c++  java
  • 3451: Tyvj1953 Normal 点分治 FFT

    国际惯例的题面:

    代价理解为重心和每个点这个点对的代价。根据期望的线性性,我们枚举每个点,计算会产生的ij点对的代价即可。
    那么,i到j的链上,i必须是第一个被选择的点。
    对于i来说,就是1/dis(i,j)。
    所以答案就是sigma(i,j) 1/(dis(i,j)+1)。
    然而这样计算是n^2的,考虑优化。
    如果我们能计算出边长为某个数值的边的数量的话,是不是就能计算答案呢?
    统计路径的题,一眼点分治。
    考虑怎样计算,我们能dfs出每个子树中距离分治重心为x的点有多少个,然后我们枚举两个点让他们取去组成路径即可。
    这显然是个卷积,FFT优化。我们补集转化,先计算全部方案,再减去本身对本身(两个点来自相同子树)的方案。
    为什么这样算复杂度正确?因为当当前分治层数一定时,所有子树的最深点的深度总和是O(n)的,并且那个log还会更小。这样分析的话发现复杂度是O(nlog^2n)。
    正常的二元关系计算方式是前缀和和当前的卷积贡献,为什么这次不能这样呢?
    给你一棵扫把形的树,一半的点形成一条链,显然你会选择扫把的重心(一边是一堆叶子,一边是链)当做重心。
    然后你发现链的那边长度为n/2,如果你对每个叶子都和链做一次卷积的话,恭喜你卡成n^2logn,不如暴力......

    代码:

     1 #include<cstdio>
     2 #include<cstring>
     3 #include<algorithm>
     4 #include<cmath>
     5 const int maxn=262145;
     6 const int inf=0x3f3f3f3f;
     7 const double pi = acos(-1.0);
     8 
     9 int tim[maxn];
    10 
    11 namespace FFT {
    12     struct Complex {
    13         double r,i;
    14         friend Complex operator + (const Complex &a,const Complex &b) { return (Complex){a.r+b.r,a.i+b.i}; }
    15         friend Complex operator - (const Complex &a,const Complex &b) { return (Complex){a.r-b.r,a.i-b.i}; }
    16         friend Complex operator * (const Complex &a,const Complex &b) { return (Complex){a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r}; }
    17     }cp[maxn];
    18     inline void FFT(Complex* dst,int n,int tpe) {
    19         for(int i=0,j=0;i<n;i++) {
    20             if( i < j ) std::swap(dst[i],dst[j]);
    21             for(int t=n>>1;(j^=t)<t;t>>=1) ;
    22         }
    23         for(int len=2;len<=n;len<<=1) {
    24             const int h = len >> 1;
    25             const Complex per = (Complex){cos(pi*tpe/h),sin(pi*tpe/h)};
    26             for(int st=0;st<n;st+=len) {
    27                 Complex w = (Complex){1.0,0.0};
    28                 for(int pos=0;pos<h;pos++) {
    29                     const Complex u = dst[st+pos] , v = dst[st+pos+h] * w;
    30                     dst[st+pos] = u + v , dst[st+pos+h] = u - v , w = w * per;
    31                 }
    32             }
    33         }
    34         if( !~tpe ) for(int i=0;i<n;i++) dst[i].r /= n;
    35     }
    36     inline void mul(int* dst,int n) {
    37         int len = 1;
    38         while( len <= ( n << 1 ) ) len <<= 1;
    39         for(int i=0;i<len;i++) cp[i] = (Complex){(double)dst[i],0.0};
    40         FFT(cp,len,1);
    41         for(int i=0;i<len;i++) cp[i] = cp[i] * cp[i];
    42         FFT(cp,len,-1);
    43         for(int i=0;i<len;i++) dst[i] = (int)(cp[i].r+0.5);
    44     }
    45 }
    46 
    47 namespace Tree {
    48     int s[maxn],t[maxn<<1],nxt[maxn<<1];
    49     int siz[maxn],mxs[maxn],ban[maxn];
    50     int su[maxn],tp[maxn];
    51     
    52     inline void addedge(int from,int to) {
    53         static int cnt = 0;
    54         t[++cnt] = to , nxt[cnt] = s[from] , s[from] = cnt;
    55     }
    56     inline void findroot(int pos,int fa,const int &fs,int &rt) {
    57         siz[pos] = 1 , mxs[pos] = 0;
    58         for(int at=s[pos];at;at=nxt[at]) if( t[at] != fa && !ban[t[at]] ) findroot(t[at],pos,fs,rt) , siz[pos] += siz[t[at]] , mxs[pos] = std::max( mxs[pos] , siz[t[at]] );
    59         if( ( mxs[pos] = std::max( mxs[pos] , fs - siz[pos]) ) <= mxs[rt] ) rt = pos;
    60     }
    61     inline void dfs(int pos,int fa,int dep,int &mxd) {
    62         mxd = std::max( mxd , dep ) , ++tp[dep];
    63         for(int at=s[pos];at;at=nxt[at]) if( t[at] != fa && !ban[t[at]] ) dfs(t[at],pos,dep+1,mxd);
    64     }
    65     inline void solve(int pos,int fs) {
    66         int root = 0 , mxd = 0 , ths ;
    67         *mxs = inf, findroot(pos,-1,fs,root) , ban[root] = 1;
    68         for(int at=s[root];at;at=nxt[at]) if( !ban[t[at]]) {
    69             ths = 0 , dfs(t[at],root,1,ths) , mxd = std::max( mxd , ths );
    70             for(int i=1;i<=ths;i++) su[i] += tp[i];
    71             FFT::mul(tp,ths);
    72             for(int i=1;i<=ths<<1;i++) tim[i] -= tp[i];
    73             memset(tp,0,sizeof(int)*(ths<<1|1));
    74         }
    75         ++*su , FFT::mul(su,mxd);
    76         for(int i=1;i<=mxd<<1;i++) tim[i] += su[i];
    77         memset(su,0,sizeof(int)*(mxd<<1|1));
    78         for(int at=s[root];at;at=nxt[at]) if( !ban[t[at]] ) solve(t[at],siz[t[at]]<siz[root]?siz[t[at]]:fs-siz[root]);
    79     }
    80 }
    81 
    82 int main() {
    83     static int n;
    84     static long double ans;
    85     scanf("%d",&n);
    86     for(int i=1,a,b;i<n;i++) scanf("%d%d",&a,&b) , ++a , ++b , Tree::addedge(a,b) , Tree::addedge(b,a);
    87     Tree::solve(1,n) , ans = n;
    88     for(int i=1;i<=n<<1;i++) ans += (long double) tim[i] / ( i + 1 );
    89     printf("%0.4Lf
    ",ans);
    90     return 0;
    91 }
    View Code


    ここでこのまま
    即使在这里就这样
    僕が消えてしまっても 誰も知らずに
    我消失不见了 谁也不会知道吧
    明日が来るのだろう
    明天依然会来临吧
    わずか 世界のひとかけらに過ぎない
    我仅仅是 这个世界的微小碎屑
    ひとりを夜が包む
    夜晚怀抱孤独的身影

  • 相关阅读:
    hibernate框架的搭建与简单实现增删改
    $.ajax();详解
    利用json实现数据传输
    利用ajax实现数据传输
    错误:Value '0000-00-00 00:00:00' can not be represented as java.sql.Timestamp;的解决
    简单使用jstl实现敏感字替换
    实用jstl实现未登录时不能绕过登录界面的效果
    简单实用jstl实现“登录|注册”
    简单实用jstl实现代码编写
    简单使用EL进行标签的替换
  • 原文地址:https://www.cnblogs.com/Cmd2001/p/8969717.html
Copyright © 2011-2022 走看看