题解 P5403 【[CTS2019]田野】

· · 题解

题目大意：给定平面上一堆线段，求一些凸包完整覆盖这些线段，求凸包周长和的最小值。

一道绝好的乱搞题。

场上对着这个题一脸懵逼，打了40分之后觉得不能就这么弃疗，就打了两个贪心拼在一起，当时想能多拿10分就算撞大运了，出来结果一看65真是激动得不行。

后来才知道我的乱搞还不算厉害，有位大神直接随机化+贪心拿到了95分……在此深受我一拜。

闲话少说进入正题。

一开始我们有n个凸包，我们的任务就是恰当地合并凸包。

我们定义某次合并的价值为：新凸包的周长-原凸包的周长和。显然这个值越小越优。

自然地产生几个问题：

1、如果某次合并的价值是负数，那这次合并是否一定要进行？

2、如果某次合并的价值是正数，那这次合并是否一定不进行？

一次合并多个凸包过于复杂，我们先从2个开始研究。

考虑这样一种贪心策略：每次选择2个凸包，如果合并的价值为负就合并，直到不存在合并的价值为负为止。

第一个问题是，在此过程中的每次合并是否一定最优？

可以证明是这样的。不过由于较（wo）为（bu）繁（hui）琐（zheng），在此略去证明。

但是这样跑出的答案就是最终答案吗？我们发现过不了最后一个样例。

我们惊奇地发现，可能出现单独合并任意两个凸包的价值都为正，但是一次性合并若干集合的价值为负的情况！

此时有一种玄学乱搞：每次随机两个凸包合并（即使价值为正），并用新集合再进行贪心合并，到不能贪心合并时继续随机，并在过程中随时更新最优答案。

多随机几次，恭喜你获得了95分的好成绩！

其实这主要是因为强数据很难造。在绝大多数测试点中，不优的合并都只需要进行1~2次，多次随机很容易就能随机出来。然而唯一没过的点需要的不优合并很多，我用我的随机程序需要跑好几分钟才能随到最优解。

所以有什么靠谱做法吗？

同样不加证明（还是因为我不会证）地给出一个引理：如果某次合并的价值为负，而从中选出任何一个子集来单独合并的价值均不为负，则这次合并一定会进行。更强的结论是，如果某次合并的价值为负，且任意子集合并的价值都大于当前集合，且不存在两个不交的子集的合并价值都为负，那么这次合并一定会进行。

注意这并不是现场讲题时ppt上给出的结论（如果任意子集合并的价值都大于当前集合就一定合并），后者可以很容易地构造数据hack掉。

因此现场讲题时给出的n^5和n^4的“每次找出全局价值最小的集合进行合并”的做法是错的！实测会在第10个点wa掉。

没错我就是一开始先打算写写n^4做法，发现各种花式wa之后强行手玩数据发现了这一点……

但是有意思的是，std的n^3做法却是正确的。实际上，std的做法就是依赖于上文给出的结论。

考虑一个对凸包的dp。

将所有点按照x第一关键字，y第二关键字排序，并从右上到左下枚举一个点。我们尝试选一个以这个点为左下角，且价值最小的凸包。

设f(i,j,k)为以点i为起点，所选的最后一条边为j~k时的最优答案。每次枚举下一个点进行转移。

但是具体怎么计算价值呢？或者换句话说，我们怎么才能知道我们新的凸包会包住哪些旧凸包呢？

有一个经典的套路叫“射线法”。对于每个凸包，我们选一个点（一般是最左下的点），向某个方向做一条射线。在dp时按逆时针顺序进行，则每选一条逆时针方向穿越某条射线的边就减去对应凸包的周长，顺时针方向穿越就加上对应周长。这样转一圈回来，所有被包进去的凸包都会被算1次贡献，而外面的凸包则全抵消掉不会算贡献。

当然这里的边界特判非常多。比如，我们不能把原先已有的凸包分割开，因此选一条边时要与所有原凸包求交判合法性；再比如，如果新凸包恰好过原来某个凸包拉射线的点，怎样才能避免这个凸包的贡献不被重算/漏算？再比如，如果拉的射线恰好过其他点怎么办……这些细节不写不知道，一写吓一跳。具体实现见代码吧。

更进一步，我们有一个喜闻乐见的优化：为什么要在dp过程中时刻保持凸性呢？显然一个不凸的方案一定不优啊。因此我们就不需要记最后一条边而是直接记最后一个点，为了防止循环转移按照极角序dp即可。

可是等等……你刚才不是说这种dp价值最小的凸包的算法不对吗？

这里非常奇妙的一点是：一旦我们找到了一个价值为负的集合，就立即将其合并！然后再从这个点开始进行新的一轮dp。如果从当前点开始没有价值为负的区间，那么我们后面的dp中就不会再考虑以它为起点的情况了！

如此一来，当我们某一次合并一个集合的时候，一定不会包含一个不以当前点为左下角的价值为负的子集（因为如果有这样的子集，它早在之前的dp中就已经被合并了），而这个集合又是以当前点为左下角的价值最小的集合，根据引理，它一定会被合并。

而为什么无法合并就可以跳过一个点呢？因为接下来的合并都一定会涉及到更靠左下的点，这自然无法使得任何一个以当前点为左下角的集合价值变小。

如此dp一次的复杂度是n^2。而每次dp要么会合并集合，要么会直接跳过一个点，这两者都只有On次，总复杂度就是n^3的了。

最后说一句，如果有人会了文中结论的证明，或者找到了把结论/std叉掉的数据，请尽快联系我直接发到uoj群里吧我实在被这个题折磨够了。

由于代码是从n^4强转过来的，所以又长又丑，大家将就着看吧……

#include<bits/stdc++.h>
using namespace std;
#define gc getchar()
#define pc putchar
#define li long long
inline li read(){
    li x = 0,y = 0,c = gc;
    while(!isdigit(c)) y = c,c = gc;
    while(isdigit(c)) x = (x << 1) + (x << 3) + (c ^ '0'),c = gc;
    return y == '-' ? -x : x;
}
#define ldb long double
int cnt,tong[510];
struct node{
    li x,y;
    node(li _x = 0,li _y = 0){x = _x;y = _y;}
}p[510];
inline bool operator < (node q,node w){
    return q.x == w.x ? q.y < w.y : q.x < w.x;
}
inline bool cmp(int q,int w){
    return p[q].x == p[w].x ? p[q].y < p[w].y : p[q].x < p[w].x;
}
inline node operator + (node q,node w){
    return node(q.x + w.x,q.y + w.y);
}
inline node operator - (node q,node w){
    return node(q.x - w.x,q.y - w.y);
}
inline li operator * (node q,node w){
    return q.x * w.y - q.y * w.x;
}
inline ldb dis(node q,node w){
    return sqrtl((q.x - w.x) * (q.x - w.x) + (q.y - w.y) * (q.y - w.y));
}
ldb ans;
struct tb{
    vector<int> a;
    ldb as;
    inline int& operator [] (int x){return a[x];}
}a[260],b[260];
int n,nw[1010],tot;
ldb f[510][510],h[510][510],ds[510][510];
int g[510][510];
int jj[510][510];
bool ff[510][510];
node nwd;
inline bool cpm(int q,int w){
    if(p[q].x > nwd.x && p[w].x < nwd.x) return 1;
    if(p[q].x < nwd.x && p[w].x > nwd.x) return 0;
    if(p[q].x == nwd.x && p[q].y < nwd.y) return 0;
    if(p[w].x == nwd.x && p[w].y < nwd.y) return 1;
    return (p[q] - nwd) * (p[w] - nwd) > 0;
}
inline bool jiao(node a,node b,node c,node d){
    node p1 = a - c,p2 = d - c,p3 = b - c;
    if(1.0l * (p1 * p2) * (p3 * p2) >= 0) return 0;
    p2 = b - a;p3 = d - a;
    return 1.0l * (p1 * p2) * (p3 * p2) > 0;
}
bool shan[260];
int pm[510];
int sx[260];
inline bool cp(int q,int w){
    return q > w;
}
bool intb[510];
int st1[1010],ft1,st2[1010],ft2;
inline ldb wk(){
    int i;
    ft1 = ft2 = 0;
    sort(nw + 1,nw + tot + 1);

    for(i = 1;i <= tot;++i){
        while(ft1 > 1 && ((p[st1[ft1]] - p[st1[ft1 - 1]]) * (p[st1[ft1]] - p[nw[i]]) >= 0)) --ft1;
        st1[++ft1] = nw[i]; 
    }
    for(i = tot;i;--i){
        while(ft2 > 1 && ((p[st2[ft2]] - p[st2[ft2 - 1]]) * (p[st2[ft2]] - p[nw[i]]) >= 0)) --ft2;
        st2[++ft2] = nw[i]; 
    }
    ldb as = 0;
    for(i = 1;i < ft1;++i) as += dis(p[st1[i]],p[st1[i + 1]]);
    for(i = 1;i < ft2;++i) as += dis(p[st2[i]],p[st2[i + 1]]);
    return as;
}
int main(){
    int i,j,k,l,u,v;
    n = read();
    for(i = 1;i <= n;++i){
        a[i].a.resize(2);
        p[++cnt].x = read();p[cnt].y = read();tong[cnt] = cnt;
        p[++cnt].x = read();p[cnt].y = read();tong[cnt] = cnt;
        a[i].as = 2 * dis(p[cnt - 1],p[cnt]);
        ans += a[i].as;
    }
    if(n == 1){printf("%.10lf\n",(double)ans);return 0;}
    sort(tong + 1,tong + cnt + 1,cmp);
    for(i = 1;i <= cnt;++i) a[tong[i] + 1 >> 1][tong[i] + 1 & 1] = i;
    for(i = 1;i <= n;++i) if(a[i][0] > a[i][1]) swap(a[i][0],a[i][1]);
    memset(tong,0,sizeof(tong));
    sort(p + 1,p + cnt + 1);
    for(i = 1;i <= cnt;++i){
        int tt = 0;
        for(j = 1;j <= cnt;++j) if(i != j) jj[i][++tt] = j;
        nwd = p[i];
        sort(jj[i] + 1,jj[i] + tt + 1,cpm);
    }
    for(i = 1;i <= cnt;++i) for(j = i;j <= cnt;++j) ds[i][j] = ds[j][i] = dis(p[i],p[j]);
    bool lss = 0,fg = 1;
    while(fg){
        fg = 0;
        for(i = 1;i < n;++i) for(j = i + 1;j <= n;++j){
            tot = 0;
            for(k = 0;k < a[i].a.size();++k) nw[++tot] = a[i][k];
            for(k = 0;k < a[j].a.size();++k) nw[++tot] = a[j][k];
            ldb nxt = wk();
            if(nxt <= a[i].as + a[j].as){
                ans += nxt - a[i].as - a[j].as;
                a[i].a.clear();a[i].a.resize(ft1 + ft2 - 2);
                for(k = 1;k < ft1;++k) a[i][k - 1] = st1[k];
                for(k = 1;k < ft2;++k) a[i][k + ft1 - 2] = st2[k];
                a[i].as = nxt;swap(a[j],a[n]);--n;fg = 1;
                goto qwq;
            }
        }
        qwq:;
    }
    for(i = 1;i <= n;++i) for(j = 0;j < a[i].a.size();++j) for(k = 1;k <= n;++k) for(l = 0;l < a[k].a.size();++l) if(i != k || j != l) ff[a[i][j]][a[k][l]] = 1;
    while(1){
        memset(g,0,sizeof(g));
        memset(intb,0,sizeof(intb));
        for(i = 1;i <= n;++i) for(j = 0;j < a[i].a.size();++j) intb[a[i][j]] = 1;
        for(i = 1;i <= cnt;++i) for(j = 1;j <= cnt;++j) if(!intb[i] || !intb[j]) ff[i][j] = 0;
        for(i = 1;i <= cnt;++i) for(j = 1;j <= cnt;++j) f[i][j] = ans;
        for(i = 1;i <= n;++i) for(j = 0;j < a[i].a.size();++j) for(k = 1;k <= n;++k) for(l = 0;l < a[k].a.size();++l) if(i != k || j != l){
            int p1 = a[i][j],p2 = a[k][l];
            if(!ff[p1][p2]) continue;
            if((p[p2] - p[p1]) * ((j != a[i].a.size() - 1 ? p[a[i][j + 1]] : p[a[i][0]]) - p[p1]) < 0){
                ff[p1][p2] = 0;continue;
            } 
            if((p[p1] - p[p2]) * ((l ? p[a[k][l - 1]] : p[a[k][a[k].a.size() - 1]]) - p[p2]) > 0){
                ff[p1][p2] = 0;continue;
            }
            for(u = 1;u <= n;++u) if(u == n || !lss){
                if((i < k || (i == k && j < l)) || !ff[p2][p1]){
                    for(v = 0;v < a[u].a.size();++v){
                        if(jiao(p[p1],p[p2],p[a[u][v]],v ? p[a[u][v - 1]] : p[a[u][a[u].a.size() - 1]])){
                            ff[p1][p2] = ff[p2][p1] = 0;break;
                        }
                    } 
                }
                if(!ff[p1][p2]) break;
                if(p1 == a[u][0]){
                    if((p[p2] - p[p1]) * node(1,-1234567890 - p[p1].y) < 0) h[p1][p2] -= a[u].as;
                } 
                else{
                    nwd = node(p[a[u][0]].x + 1,-1234567890);
                    li q1 = (p[p1] - p[a[u][0]]) * node(1,-1234567890 - p[a[u][0]].y),q2 = (p[p2] - p[a[u][0]]) * node(1,-1234567890 - p[a[u][0]].y);
                    if(q1 > 0 && q2 < 0){
                        if(1.0l * ((nwd - p[p1]) * (p[p2] - p[p1])) * ((p[a[u][0]] - p[p1]) * (p[p2] - p[p1])) < 0) h[p1][p2] -= a[u].as;
                    } 
                    else if(q1 < 0 && q2 > 0){
                        if(1.0l * ((nwd - p[p1]) * (p[p2] - p[p1])) * ((p[a[u][0]] - p[p1]) * (p[p2] - p[p1])) < 0) h[p1][p2] += a[u].as;   
                    } 
                }
            } 
        }
        ldb nwa = ans;int wz = 0;
        for(i = 1;i <= n;++i) sx[i] = a[i][0];
        sort(sx + 1,sx + n + 1,cp);
        for(i = 1;i <= n;++i){
            int p1 = sx[i],p2;
            if(lss && p[a[n][0]] < p[p1]) continue; 
            pm[p1] = cnt;
            for(j = 1;j < cnt;++j) pm[jj[p1][j]] = j;
            for(j = 1;j < cnt && p[jj[p1][j]].x >= p[p1].x;++j){
                p2 = jj[p1][j];
                if(p[p2].x == p[p1].x && p[p2].y < p[p1].y) break;
                if(!ff[p1][p2]) continue;
                if(f[p1][p2] > ds[p1][p2] + h[p1][p2]){
                    f[p1][p2] = ds[p1][p2] + h[p1][p2];
                    g[p1][p2] = p1;
                }
            }
            for(j = 1;j < cnt && p[jj[p1][j]].x >= p[p1].x;++j){
                p2 = jj[p1][j];
                if(p[p2].x == p[p1].x && p[p2].y < p[p1].y) break;
                for(k = 1;k <= cnt;++k) if(ff[p2][k] && pm[k] > pm[p2]){
                    if(f[p1][k] > f[p1][p2] + ds[p2][k] + h[p2][k]){
                        f[p1][k] = f[p1][p2] + ds[p2][k] + h[p2][k];
                        g[p1][k] = p2;
                    }
                }
            }
            if(f[p1][p1] < -1e-8){
                nwa = f[p1][p1];wz = p1;break;
            } 
        }
        if(nwa >= -1e-8) break;
        ans += nwa;
        nw[tot = 1] = wz;
        for(i = g[wz][wz];i != wz;i = g[wz][i]) nw[++tot] = i;
        nw[++tot] = wz;
        int nxtn = 0;
        for(i = 1;i <= n;++i){
            nwd = node(p[a[i][0]].x + 1,-1234567890);
            bool inn = 0;
            for(j = 1;j < tot;++j){
                if(p[a[i][0]].x == p[nw[j]].x && p[a[i][0]].y == p[nw[j]].y){
                    inn = 1;break;
                }
                if(jiao(p[nw[j]],p[nw[j + 1]],p[a[i][0]],nwd)) inn ^= 1;
            }
            if(inn){
                nwa += a[i].as;u = i;
                for(int p1 = 1;p1 <= cnt;++p1) for(int p2 = 1;p2 <= cnt;++p2) if(ff[p1][p2]){
                    if(p1 == a[u][0]){
                        if((p[p2] - p[p1]) * node(1,-1234567890 - p[p1].y) < 0) h[p1][p2] += a[u].as;
                    } 
                    else{
                        nwd = node(p[a[u][0]].x + 1,-1234567890);
                        li q1 = (p[p1] - p[a[u][0]]) * node(1,-1234567890 - p[a[u][0]].y),q2 = (p[p2] - p[a[u][0]]) * node(1,-1234567890 - p[a[u][0]].y);
                        if(q1 > 0 && q2 < 0){
                            if(1.0l * ((nwd - p[p1]) * (p[p2] - p[p1])) * ((p[a[u][0]] - p[p1]) * (p[p2] - p[p1])) < 0) h[p1][p2] += a[u].as;
                        } 
                        else if(q1 < 0 && q2 > 0){
                            if(1.0l * ((nwd - p[p1]) * (p[p2] - p[p1])) * ((p[a[u][0]] - p[p1]) * (p[p2] - p[p1])) < 0) h[p1][p2] -= a[u].as;   
                        } 
                    }
                }
            } 
            else b[++nxtn] = a[i];
        }
        for(i = 1;i <= nxtn;++i) a[i] = b[i];
        n = nxtn + 1;
        a[n].as = nwa;a[n].a.clear();a[n].a.resize(tot - 1);
        for(i = tot;i > 1;--i) a[n][tot - i] = nw[i];lss = 1;
    }
    printf("%.10lf\n",(double)ans);
    return 0;
}