BZOJ:3911: SGU383 Caravans（三角剖分）

原题链接：http://www.lydsy.com/JudgeOnline/problem.php?id=3911

直接求最小生成树显然边太多，考虑少用点边。

连出来的边肯定是没相交的，我们需要做一下三角剖分，求出可能的待选边。

这个很资瓷：http://www.geom.uiuc.edu/~samuelp/del_project.html

三角剖分大约就是分治到3个点以内直接两两连边，合并两个块的时候需要先求出下凸包，找到最底下的一条边先连起来，然后找这条边与其他未处理的点能构成的最小圆（没有其他点在这个圆内）。删掉相交的边然后不断把底边往上抬（这里需要结合图来理解）。

最后有一些细节诸如求凸包的时候三点共线要保留什么的。一开始还在纠结删边的时候点在边上怎么处理，后来才发现并不会有这种情况。

#include<cmath>
#include<cstdio>
#include<algorithm>
#define eps 1e-12
#define MN 100001
#define ll long long
#define ld double
using namespace std;

int read_p,read_ca;
inline int read(){
    read_p=0;read_ca=getchar();
    while(read_ca<'0'||read_ca>'9') read_ca=getchar();
    while(read_ca>='0'&&read_ca<='9') read_p=read_p*10+read_ca-48,read_ca=getchar();
    return read_p;
}
struct po{int x,y;};
struct na{po p;int n;}p[MN];
struct e{int x,y,ne,pre;}b[MN<<5];
struct _e{int y,ne;ll z;}_b[MN];
struct bi{int x,y;ll z;}B[MN<<5];
bool operator < (bi a,bi b){return a.z<b.z;}
bool operator < (na a,na b){return a.p.x==b.p.x?a.p.y<b.p.y:a.p.x<b.p.x;}
int n,m,num=1,l[MN],st[MN],top,ct[MN],NUM=0,x,y,_l[MN],_num=0,fa[MN],f[MN][20],de[MN];
ll v[MN][20];
inline void in(int x,int y,ll z){_b[++_num].y=y;_b[_num].z=z;_b[_num].ne=_l[x];_l[x]=_num;}
inline void _add(int x,int y){b[++num].x=x;b[num].y=y;b[num].ne=l[x];b[l[x]].pre=num;l[x]=num;}
inline void add(int x,int y){_add(x,y);_add(y,x);}
inline void _del(int i){if(l[b[i].x]==i)l[b[i].x]=b[i].ne;else{b[b[i].pre].ne=b[i].ne;if(b[i].ne)b[b[i].ne].pre=b[i].pre;}}
inline void del(int i){_del(i);_del(i^1);}
inline int max(int a,int b){return a>b?a:b;}
ll cha(po a,po b,po c){return 1LL*(b.x-a.x)*(c.y-b.y)-1LL*(c.x-b.x)*(b.y-a.y);}
bool _cha(po a,po b,po c){return cha(a,b,c)<=0;}
inline bool ju(po a,po b,po c,po d){return (_cha(a,b,c)^_cha(a,b,d))&&(_cha(c,d,a)^_cha(c,d,b));}
ll sqr(int &x){return 1LL*x*x;}
ld sqr(ld x){return x*x;}
bool incir(po a,po b,po c,po d){
    ll m=(1LL*a.x*(b.y-c.y)+1LL*b.x*(c.y-a.y)+1LL*c.x*(a.y-b.y))<<1;
    ld x=(ld)(sqr(a.x)*(b.y-c.y)+sqr(b.x)*(c.y-a.y)+sqr(c.x)*(a.y-b.y)-1LL*(a.y-b.y)*(b.y-c.y)*(c.y-a.y))/m;
    ld y=(ld)(1LL*(a.x-b.x)*(b.x-c.x)*(c.x-a.x)-sqr(a.y)*(b.x-c.x)-sqr(b.y)*(c.x-a.x)-sqr(c.y)*(a.x-b.x))/m;
    return sqrt(sqr(x-a.x)+sqr(y-a.y))+eps>sqrt(sqr(x-d.x)+sqr(y-d.y));
}
void work(int L,int R){
    if (R-L<=2){
        for (int i=L;i<=R;i++) ct[p[i].n]=i;
        for (int i=L;i<R;i++)
        for (int j=i+1;j<=R;j++) add(i,j);
        return;
    }
    int mid=L+R>>1;
    work(L,mid);work(mid+1,R);
    top=0;
    for (int i=L;i<=R;st[++top]=i++)
    while (top>1&&cha(p[st[top-1]].p,p[st[top]].p,p[i].p)<0) top--;
    int Ld,Rd,id;
    for (int i=1;i<top;i++)
    if (st[i]<=mid&&st[i+1]>mid) Ld=st[i],Rd=st[i+1];
    add(Ld,Rd);
    for(;;){
        id=-1;
        for (int i=l[Ld];i;i=b[i].ne)
        if (b[i].y!=Rd)
        if (cha(p[b[i].y].p,p[Ld].p,p[Rd].p)>0&&(id==-1||incir(p[Ld].p,p[Rd].p,p[id].p,p[b[i].y].p))) id=b[i].y;
        for (int i=l[Rd];i;i=b[i].ne)
        if (b[i].y!=Ld)
        if (cha(p[Ld].p,p[Rd].p,p[b[i].y].p)>0&&(id==-1||incir(p[Ld].p,p[Rd].p,p[id].p,p[b[i].y].p))) id=b[i].y;
        if (id==-1) break;
        
        if (id<=mid){
            for (int i=l[Ld];i;i=b[i].ne)
            if (b[i].y!=Rd)
            if (ju(p[id].p,p[Rd].p,p[b[i].y].p,p[Ld].p)) del(i);
            Ld=id;
        }else{
            for (int i=l[Rd];i;i=b[i].ne)
            if (b[i].y!=Ld)
            if (ju(p[id].p,p[Ld].p,p[b[i].y].p,p[Rd].p)) del(i);
            Rd=id;
        }
        add(Ld,Rd);
    }
}
int gf(int x){return fa[x]==x?x:fa[x]=gf(fa[x]);}
void dfs(int x,int _f){
    for (int i=1;i<20;i++) f[x][i]=f[f[x][i-1]][i-1],v[x][i]=max(v[x][i-1],v[f[x][i-1]][i-1]);
    for (int i=_l[x];i;i=_b[i].ne)
    if (_b[i].y!=_f) de[_b[i].y]=de[x]+1,f[_b[i].y][0]=x,v[_b[i].y][0]=_b[i].z,dfs(_b[i].y,x);
}
ll Mavis(int x,int y){
    ll mmh=0;
    if (de[x]<de[y]) swap(x,y);
    for (int i=19;i>=0;i--)
    if (de[y]<=de[f[x][i]]){
        if (v[x][i]>mmh) mmh=v[x][i];
        x=f[x][i];
    }
    if (x==y) return mmh;
    for (int i=19;i>=0;i--)
    if (f[y][i]!=f[x][i]){
        if (v[x][i]>mmh) mmh=v[x][i];
        if (v[y][i]>mmh) mmh=v[y][i];
        x=f[x][i];y=f[y][i];
    }
    if (v[x][0]>mmh) mmh=v[x][0];
    if (v[y][0]>mmh) mmh=v[y][0];
    return mmh;
}
int main(){
    register int i,j;
    n=read();
    for (i=0;i<n;i++) p[i].p.x=read(),p[i].p.y=read(),p[i].n=i+1,fa[i]=i;
    sort(p,p+n);
    work(0,n-1);
    for (i=0;i<n;i++)
    for (j=l[i];j;j=b[j].ne) if (b[j].y>i) B[NUM].x=i,B[NUM].y=b[j].y,B[NUM++].z=sqr(p[i].p.x-p[b[j].y].p.x)+sqr(p[i].p.y-p[b[j].y].p.y);
    sort(B,B+NUM);
    for (i=0;i<NUM;i++){
        if (p[B[i].x].n>p[B[i].y].n) swap(B[i].x,B[i].y);
        int x=gf(B[i].x),y=gf(B[i].y);
        if (x==y) continue;
        fa[x]=y;
        in(B[i].x,B[i].y,B[i].z);in(B[i].y,B[i].x,B[i].z);
    }
    dfs(0,0);
    m=read();
    while (m--){
        x=ct[read()];y=ct[read()];
        printf("%.6lf
",sqrt(Mavis(x,y)));
    }
}

---

upd 2017/12/26

　　代码被hack了，当初应该是对着n+e的blog学的吧，所以他n^2我也n^2了……

　　最近已经无力捣鼓优化惹，直接贴一下别人blog吧：https://www.cnblogs.com/darklove/p/8056097.html