半平面交对偶转凸包树分治归并排序求闵可夫斯基和用李超线段树维护答案。
题意:给一颗树,每条边的边权是一个关于t的一次函数,问你在每个t下树的直径是多少。
如果t是定值,那么就是求直径。求直径的方法除了dfs两遍,还有一个方法是树分治,类似分治法求最大子段和。
考虑树分治,找到重心,枚举所有从重心延伸到每个点的一次函数y=ax+b。这个是比较套路的,先统计每个子树的信息嘛。
统计信息的时候我们先做两个优化:
1、首先能够想到因为题目给出的a,b均正,又让求最大直径。所以只有从重心延伸到叶子节点的一次函数才是有用的。此处可以进行剪枝,仅保留到叶子节点时的信息。
2、每个叶节点提供的是一个at+b的一次函数,随着t的变化,有些曲线不可能成为最优曲线,所以此处可以先做个半平面交。就是仅保留从上方无限远往下看可以看到的曲线。
肯定不能
合并所有子树上的一次函数,即使有着两个小优化,也没有根本性的解决问题,还是会被菊花树卡死。
考虑怎么在合并上做优化。
现在的问题就是说每个重心有若干个子树,每个子树上都储存了一个半平面交构成的凸包,如何合并信息?
先假设对于每个分治重心,只有两个子树的情况。
黑科技:闵可夫斯基和
定义平面上的点的加法,若有两点
,
则
定义平面上的图形A,B为A,B图形内所有的点组成的集合。
则可以定义两个图形在平面内的加法
这个东西虽然看着不太像能够快速求出来的样子,不过如果图形A和图形B是两个凸多边形,就有了一个线性的求法。
定理:两个凸多边形的闵可夫斯基和仍是一个凸多边形。
定理:两个凸多边形的闵可夫斯基和上顶点的个数总等于这两个凸多边形每条有向边向量角的种类数。
定理:两个凸多边形的闵可夫斯基和的周长总等于这两个凸多边形的周长之和。
粉色图形是两个蓝色凸多边形的闵可夫斯基和。
首先把这两个凸多边形按照顺时针或者逆时针的顺序整理好,取两个图形的第一个点作为闵可夫斯基和的起点,按级角序依次弹出每一条边。
由于凸包平时就已经是按照级角序排好储存的。所以求两个凸包的闵可夫斯基和就类似一次归并排序。
typedef vector<point> pointset;
typedef vector<point> convex;
void convex_minkowski_sum(const convex &A,const convex &B,convex & ans)
{
point pos=A[0]+B[0];
ans.clear();
ans.push_back(pos);
int p1=1,p2=1;
while(p1<=A.size()&&p2<=B.size())
{
if(dcmp(cross(A[p1%A.size()]-A[p1-1],B[p2%B.size()]-B[p2-1]))>=0)
{
pos=pos+A[p1%A.size()]-A[p1-1];
++p1;
}
else
{
pos=pos+B[p2%B.size()]-B[p2-1];
++p2;
}
ans.push_back(pos);
}
while(p1<=A.size())
{
pos=pos+A[p1%A.size()]-A[p1-1];
++p1;
ans.push_back(pos);
}
while(p2<=B.size())
{
pos=pos+B[p2%B.size()]-B[p2-1];
++p2;
ans.push_back(pos);
}
ans.pop_back();
}
有了闵可夫斯基和,我们能够快速合并两个子树之间的所有一次函数。
只要对这两个子树半平面交构成的凸包做一次闵可夫斯基和,那么显然这两个子树上所有一次函数两两组合的最优曲线也一定在和凸包上。
然后你发现,计算闵可夫斯基和需要的是点围成的凸包,而半平面交是线组成的,虽然能够先求的所有的交点扔进去跑,但是一来一回的转化不够优美,精度也飞到不知道哪去了。
所以用到了第二个黑科技。
黑科技:半平面交对偶转凸包
其实就是做一种对偶变换。
对于直线y=kx+b,可以将其映射到平面上的一点
,点p就是直线的对偶点。
那么对于若干条直线组成的上凸壳,等于它们对偶点组成的下凸壳。
对于若干条直线组成的下凸壳,等于它们对偶点组成的上凸壳。
这样使用起来容易出错,因为对偶之后凸壳的上下是颠倒的。实际使用中往往采用
,作为直线y=kx+b的对偶点。
这样对偶凸壳的上下关系也跟原凸壳对应了起来,不容易出错。
一开始求半平面交的时候就转化为对偶点组成的“点”凸包,就避开了后面求闵可夫斯基和时“点”凸包和“边”凸包的相互转化,可以减小精度误差。
有了闵可夫斯基和跟半平面交对偶转凸包,可以很轻松的合并两个子树组成的凸壳。
至此,总结一下我们要干的事:
1、树分治求出每个子树延伸到叶子节点中所有的一次函数y=kx+b,将其对偶化求凸包。
2、假设每次分治只用合并两个子树,那么求出这两个子树的闵可夫斯基和,和凸包上的每个坐标点p(k,-b),都对应一条曲线y=kx+b,这条曲线可能出现在最优凸包上。
3、将所有可能出现的曲线都拿出来再做一个大凸包,顺着for过去依次查询每个时刻即可。
到这里已经可以使用官方题解中的思路,将多叉树转化为二叉树,然后使用树的边分治来解决此题
但是我偏要用点分治!
点分治比起边分治最大的问题在于,边分治是两两合并,而点分治面临一堆东西之间的相互合并,如果处理不好,还是
合并。
闵可夫斯基和类似单次归并排序,归并排序它能干嘛呢?
归并排序能求逆序对!!!
归并排序是一个天然的
的优化,它是如何只用nlog的代价枚举求出两两之间的逆序关系的呢?
我们画出归并树,发现任意两个元素的逆序信息将在这两个元素在归并树上的LCA节点被处理。
我们将这若干个子树提供的凸壳看成是元素。
仿照归并排序求逆序对,我们写一个归并排序求闵可夫斯基和,那么任意两个子树凸壳的合并信息将会在处理这两个子树元素的归并树上的LCA节点被更新到。
这样就做到了只用
的代价,就处理了n个元素之间两两的合并关系。
至于最后画这个大凸壳,理论上讲没有必要用李超树,使用单调栈维护即可。但是还要重新排序,而李超树在处理直线是
复杂度,(李超树处理线段时,复杂度会退化为
)
反正排序也带log,既然使用单调栈也还得带log,我就不麻烦了直接上李超树。
#include<bits/stdc++.h>
using namespace std;
long double zero;
struct Lin
{
long double k,b;
long double val(long double x)
{
return k*x+b;
}
Lin(long double _k=0,long double _b=0)
{
k=_k;
b=_b;
}
};
const long double EPS=1e-7;
const long double INF=1e18;
struct LiChao_Segmenttree
{
int dcmp(long double x)
{
return x>EPS?1:x<EPS?-1:0;
}
long double cross_x(const Lin&A,const Lin&B)
{
return (B.b-A.b)/(A.k-B.k);
}
struct tree_node
{
int l,r;
bool vis;
Lin lin;
} t[4*1000005];
void build(int l,int r,int root=1)
{
t[root].l=l;
t[root].r=r;
t[root].vis=false;
if(l!=r)
{
int mid=(l+r)>>1;
int ch=root<<1;
build(l,mid,ch);
build(mid+1,r,ch+1);
}
}
void insert(Lin x,int root=1)
{
if(!t[root].vis)
{
t[root].vis=true;
t[root].lin=x;
return;
}
if(dcmp(t[root].lin.val(t[root].l)-x.val(t[root].l))>=0&&dcmp(t[root].lin.val(t[root].r)-x.val(t[root].r))>=0)
{
return;
}
if(dcmp(t[root].lin.val(t[root].l)-x.val(t[root].l))<0&&dcmp(t[root].lin.val(t[root].r)-x.val(t[root].r))<0)
{
t[root].lin=x;
return;
}
int mid=(t[root].l+t[root].r)>>1;
int ch=root<<1;
if(dcmp(cross_x(x,t[root].lin)-mid)<=0)
{
if(x.k<t[root].lin.k)
{
insert(x,ch);
}
else
{
insert(t[root].lin,ch);
t[root].lin=x;
}
}
else
{
if(x.k>t[root].lin.k)
{
insert(x,ch+1);
}
else
{
insert(t[root].lin,ch+1);
t[root].lin=x;
}
}
}
void clear(int root=1)
{
t[root].vis=false;
if(t[root].l!=t[root].r)
{
int ch=root<<1;
if(t[ch].vis)
{
clear(ch);
}
if(t[ch+1].vis)
{
clear(ch+1);
}
}
return;
}
long double get_val(int x,int root=1)
{
if(!t[root].vis)return -INF;
long double ret=t[root].lin.val(x);
if(t[root].l==t[root].r)return ret;
int ch=root<<1;
int mid=(t[root].l+t[root].r)>>1;
if(x<=mid)
{
return max(ret,get_val(x,ch));
}
else
{
return max(ret,get_val(x,ch+1));
}
}
} seg;
const int MAXN=100005,MAXM=100005;
const long double PI = 4.0 * atan(1.0);
struct point
{
long double x,y;
point (long double _x = 0,long double _y = 0)
{
x = _x;
y = _y;
}
};
typedef point vect;
int dcmp(long double x)
{
return x<-EPS?-1:(x>EPS?1:0);
}
bool operator < (const point& a,const point& b)
{
if(a.x!=b.x)
return a.x<b.x;
return a.y<b.y;
}
bool operator == (const point& a,const point& b)
{
return !dcmp(a.x-b.x)&&!dcmp(a.y-b.y);
}
bool operator != (const point& a,const point& b)
{
return dcmp(a.x-b.x)||dcmp(a.y-b.y);
}
vect operator + (vect a,vect b)
{
return vect (a.x+b.x,a.y+b.y);
}
vect operator - (point a,point b)
{
return vect (a.x-b.x,a.y-b.y);
}
vect operator * (vect a,long double b)
{
return vect(a.x*b,a.y*b);
}
long double cross(vect a,vect b)
{
return a.x*b.y-b.x*a.y;
}
typedef vector<point> pointset;
typedef vector<point> convex;
void convex_hull_points(pointset &p,convex &ans)
{
sort(p.begin(),p.end());
ans.clear();
for(int i=0; i<p.size(); ++i)
{
while(ans.size()>1&&cross(ans[ans.size()-1]-ans[ans.size()-2],p[i]-ans[ans.size()-2])<=0)ans.pop_back();
ans.push_back(p[i]);
}
return;
}
pointset ch_pointset[MAXN];
convex ch_convex[MAXN];
int n,m,g[MAXN],tot,gravity_size[MAXN],gravity_max[MAXN],gravity,siz,u,v,c,du[MAXN],ch_tot;
long double k,b;
bool div_vis[MAXN];
struct edges
{
int to;
long double v1;
long double v2;
int next;
} e[MAXM<<1];
void add_edge(int from,int to,long double v1,long double v2)
{
e[++tot].to=to;
e[tot].v1=v1;
e[tot].v2=v2;
e[tot].next=g[from];
g[from]=tot;
return;
}
void get_gravity(int x,int fa)
{
gravity_size[x]=1;
gravity_max[x]=-1;
for(int i=g[x]; i; i=e[i].next)
{
if(e[i].to!=fa&&!div_vis[e[i].to])
{
get_gravity(e[i].to,x);
if(gravity_max[x]==-1||gravity_max[x]<gravity_size[e[i].to])
{
gravity_max[x]=gravity_size[e[i].to];
}
gravity_size[x]+=gravity_size[e[i].to];
}
}
gravity_max[x]=max(gravity_max[x],siz-gravity_size[x]);
if(gravity==-1||gravity_max[gravity]>gravity_max[x])gravity=x;
}
void get_size(int x,int fa)
{
++siz;
for(int i=g[x]; i; i=e[i].next)
{
if(e[i].to!=fa&&!div_vis[e[i].to])
{
get_size(e[i].to,x);
}
}
}
void update(int x,int fa,long double nowk,long double nowb,pointset &Data)
{
if(du[x]==1)
{
Data.push_back(point(nowk,nowb));
}
for(int i=g[x]; i; i=e[i].next)
{
if(e[i].to!=fa&&!div_vis[e[i].to])
{
update(e[i].to,x,nowk+e[i].v1,nowb+e[i].v2,Data);
}
}
return;
}
convex merge_sort(int l,int r)
{
convex ret;
if(l>r)return ret;
if(l==r)return ch_convex[l];
int mid=(l+r)>>1;
convex L=merge_sort(l,mid);
convex R=merge_sort(mid+1,r);
if(L.size()==0)return R;
if(R.size()==0)return L;
convex sum;
pointset ps;
point pos=L[0]+R[0];
sum.push_back(pos);
int p1=1,p2=1;
while(p1<=L.size()&&p2<=R.size())
{
if(dcmp(cross(L[p1%L.size()]-L[p1-1],R[p2%R.size()]-R[p2-1]))>=0)
{
pos=pos+L[p1%L.size()]-L[p1-1];
p1++;
}
else
{
pos=pos+R[p2%R.size()]-R[p2-1];
p2++;
}
sum.push_back(pos);
}
while(p1<=L.size())
{
pos=pos+L[p1%L.size()]-L[p1-1];
p1++;
sum.push_back(pos);
}
while(p2<=R.size())
{
pos=pos+R[p2%R.size()]-R[p2-1];
p2++;
sum.push_back(pos);
}
for(auto &i:sum)
{
seg.insert(Lin(i.x,-i.y));
zero=max(zero,-i.y);
}
p1=0,p2=0;
while(p1<L.size()&&p2<R.size())
{
if(L[p1]<R[p2])
{
while(ret.size()>1&&cross(ret[ret.size()-1]-ret[ret.size()-2],L[p1]-ret[ret.size()-2])<=0)ret.pop_back();
ret.push_back(L[p1++]);
}
else
{
while(ret.size()>1&&cross(ret[ret.size()-1]-ret[ret.size()-2],R[p2]-ret[ret.size()-2])<=0)ret.pop_back();
ret.push_back(R[p2++]);
}
}
while(p1<L.size())
{
while(ret.size()>1&&cross(ret[ret.size()-1]-ret[ret.size()-2],L[p1]-ret[ret.size()-2])<=0)ret.pop_back();
ret.push_back(L[p1++]);
}
while(p2<R.size())
{
while(ret.size()>1&&cross(ret[ret.size()-1]-ret[ret.size()-2],R[p2]-ret[ret.size()-2])<=0)ret.pop_back();
ret.push_back(R[p2++]);
}
return ret;
}
void tree_div(int x)
{
siz=0;
get_size(x,-1);
gravity=-1;
get_gravity(x,-1);
div_vis[gravity]=true;
ch_tot=0;
for(int i=g[gravity]; i; i=e[i].next)
{
if(!div_vis[e[i].to])
{
ch_pointset[++ch_tot].clear();
update(e[i].to,gravity,e[i].v1,e[i].v2,ch_pointset[ch_tot]);
}
}
for(int i=1; i<=ch_tot; ++i)
{
convex_hull_points(ch_pointset[i],ch_convex[i]);
}
merge_sort(1,ch_tot);
for(int i=g[gravity]; i; i=e[i].next)
{
if(!div_vis[e[i].to])
{
tree_div(e[i].to);
}
}
return;
}
int main()
{
scanf("%d %d",&n,&m);
if(n==1)
{
for(int i=1; i<=m; ++i)
{
printf("0%c",i==m?'\n':' ');
}
return 0;
}
if(n==2)
{
scanf("%d %d %Lf %Lf",&u,&v,&k,&b);
for(int i=1; i<=m; ++i)
{
printf("%.0Lf%c",k*(i-1)+b,i==m?'\n':' ');
}
return 0;
}
seg.build(1,m);
for(int i=1; i<n; ++i)
{
scanf("%d %d %Lf %Lf",&u,&v,&k,&b);
b=-b;
add_edge(u,v,k,b);
add_edge(v,u,k,b);
du[u]++;
du[v]++;
}
tree_div(1);
printf("%.0Lf",zero);
for(int i=1; i<m; ++i)
{
printf(" %.0Lf",max(seg.get_val(i),(long double)0),i==m?'\n':' ');
}
printf("\n");
return 0;
}
代码长主要是都是板子拼在一起...实际并没有手写什么东西
京公网安备 11010502036488号