UVALive 4986 Dome of Circus(三分、凸包、凸性函数)

题目链接；
UVALive 4986 Dome of Circus
题意：

在空间中给n个点的坐标，并且每个点都是在z轴的正半平面内，以原点(0,0,0)为底面圆圆心，画一个圆锥，要求圆锥要包含所有的点，点可以在圆锥面上，求最小的圆锥体积的底面圆半径和高。
数据范围：1≤n≤10000,并且所有点的坐标的绝对值都不超过10000
分析；
首先把三维的点转成二维的，因为是圆锥覆盖，所以我们只需要知道点到圆锥轴的距离和点的纵坐标就好了。那么可得：

(x,y,z)→(x2+y2−−−−−−√,z)

这样子一来我们就把所有的点限制在第一象限内。
我们来考虑当第一象限内的一个点(x,y)恰好在圆锥面上时的各个变量之间的关系，设此时圆锥的底面圆半径和高分别为r和h。通过相似三角形可以得到：

h−zx2+y2−−−−−−√=hr

化简一下可得:

h=z∗rr−x2+y2−−−−−−√

我们带进圆锥的体积公式可得：

v=13∗π∗h∗r2=13∗π∗z∗r3r−x2+y2−−−−−−√

我们对v求两次导可得（抛开系数13∗π∗z不看，并令t=x2+y2−−−−−−√）：

v′=r2∗(2∗r−3∗t)(r−t)2

v′′=2∗r∗(r2+3∗t2−3∗r∗t)(r−t)3=2r∗[(r−32t)2+14t2](r−t)3

于是我们可以发现v′′其实是一个恒大于0的函数，所以函数v′是单调递增的。同时我们可以发现函数v′是先负后正的。那么函数v关于r就是先增后减的，是个凸函数。

一元可微函数在某个区间上是凸的，当且仅当它的导数在该区间上单调不减。

于是我们可以三分枚举底面圆半径r，同时根据：

h=z∗rr−x2+y2−−−−−−√

对于每个点如果r确定了，那么最小体积圆锥的h也就确定了，为了满足所有点都覆盖，需要取hmax。比较r=mid和r=midmid时的体积，取最小即可。

但是我们也可以换个思路。回到最初的问题，我们是想在第一象限内画一条斜率为负的直线，并且假设和横轴、纵轴分别交于(r,0),(0,h)，如何使得h∗r2最小？
在看看上面的推导根据v的一阶导数可知当v′=0时，可以使得v最小，也就是r=32∗t,其中t=x2+y2−−−−−−√,那么也就是说对于每个点如果单独考虑的话，我们可以直接取r=32∗t,其中t=x2+y2−−−−−−√，h=z∗rr−x2+y2√使得这个圆锥体积最小，设这时直线为l。但是因为要考虑所有的点，所以选取的点肯定是凸包上的点，那么凸包上的点和这条直线l的关系有这几种：

我们通过求凸包边和横轴的交点和直线l和横轴的交点比较可以确定直线l和凸包边的关系是L1,L2,L3中的哪种。我们考虑l和p点下方的边比较。

• 如果是L1或者L3,那么最小圆锥的母线只能选择L
• 如果是L2，那么最小圆锥的母线可以选择L2，并且这种情况是最优的

我们通过求凸包，枚举凸包边，比较l和边与横轴交点，对每个顶点求能覆盖的最小圆锥，也可已解决的。
但是写起来太麻烦了。。。。而且边界情况比较多。
下面的两个代码都是AC了的。

#include <iostream>#include <cstdio>#include <cstring>#include <string>#include <algorithm>#include <climits>#include <cmath>#include <ctime>#include <cassert>#define IOS ios_base::sync_with_stdio(0); cin.tie(0);using namespace std;typedef long long ll;const int MAX_N = 100010;const double eps = 1e-8;int n;inline int sgn(double x) {    if(fabs(x) <= eps) return 0;    else if(x > 0.0) return 1;    else return -1;}struct Point {    double x, y;    Point() {}    Point(double _x, double _y): x(_x), y(_y) {}}point[MAX_N];inline double GetHeight(double r){    double res = -1.0;    for(int i = 0; i < n; ++i) {        double h = point[i].y * r / (r - point[i].x);        if(sgn(h - res) > 0) res = h;    }    return res;}int main(){    while(~scanf("%d", &n)) {        double MinR = -1.0;        for(int i = 0; i < n; ++i) {            double a, b, c;            scanf("%lf%lf%lf", &a, &b, &c);            point[i] = Point(sqrt(a * a + b * b), c);            if(sgn(point[i].x - MinR) > 0) MinR = point[i].x;        }        double low = MinR, high = 1e10, mid, midmid, v1, v2;        for(int i = 0; i < 100; ++i) {            mid = (low + high) / 2.0;            midmid = (mid + high) / 2.0;            v1 = GetHeight(mid) * mid * mid;            v2 = GetHeight(midmid) * midmid * midmid;            if(sgn(v2 - v1) < 0) low = mid;            else high = midmid;        }         printf("%.3lf %.3lf\n", GetHeight(mid), mid);    }    return 0;}

#include <iostream>#include <cstdio>#include <cstring>#include <string>#include <algorithm>#include <climits>#include <cmath>#include <ctime>#include <cassert>#include <vector>#define IOS ios_base::sync_with_stdio(0); cin.tie(0);using namespace std;typedef long long ll;const int MAX_N = 10010;const double eps = 1e-8;inline int sgn(double x) {    if(fabs(x) <= eps) return 0;    else if(x > 0) return 1;    else return -1;}struct Point {    double x, y;    Point() {}    Point(double _x, double _y): x(_x), y(_y) {}    Point operator - (const Point& rhs) const {        return Point(x - rhs.x, y - rhs.y);    }     Point operator + (const Point& rhs) const {        return Point(x + rhs.x, y + rhs.y);    }    Point operator * (const double& d) const {        return Point(x * d, y * d);    }    Point operator / (const double& d) const {        return Point(x / d, y / d);    }    double dis(const Point& rhs) const {        return sqrt((x - rhs.x) * (x - rhs.x) + (y - rhs.y) * (y - rhs.y));    }    double cross(const Point& rhs) const {        return x * rhs.y - y * rhs.x;    }    double dot(const Point& rhs) const {        return x * rhs.x + y * rhs.y;    }}point[MAX_N], vertex[MAX_N];int n;Point center = Point(0, 0);bool cmp_pola_angle(Point a, Point b){    double res = (a - center).cross(b - center);    if(sgn(res) != 0) return sgn(res) > 0;    else return a.dis(center) < b.dis(center);}int GetConvex(){    sort(point, point + n, cmp_pola_angle);    int k = 0;    for (int i = 0; i < n; ++i) {        while(k > 1 && sgn((vertex[k - 1] - vertex[k - 2]).cross(point[i] - vertex[k - 1])) <= 0) {            k--;        }        vertex[k++] = point[i];    }    int m = k;    for (int i = n - 2; i >= 0; --i) {        if(k > m && sgn((vertex[k - 1] - vertex[k - 2]).cross(point[i] - vertex[k - 1])) <= 0) {            k--;        }        vertex[k++] = point[i];    }    if(k > 1) k--;    return k;}inline double GetX(Point a, Point b){    return (a.x * b.y - a.y * b.x) / (b.y - a.y);}inline double GetY(Point a, Point b){    return (a.y * b.x - a.x * b.y) / (b.x - a.x);}int main(){    while(~scanf("%d", &n)) {        for(int i = 0; i < n; ++i) {            double a, b, c;            scanf("%lf%lf%lf", &a, &b, &c);            point[i] = Point(sqrt(a * a + b * b), c);        }        int m = GetConvex();        double r, h, v;        if(m == 1) {            printf("%.3lf %.3lf\n", 3 * vertex[0].y, 1.5 * vertex[0].x);            continue;        } else if(m == 2) {            if(vertex[0].y > vertex[1].y) swap(vertex[0], vertex[1]);            double tmp = GetX(vertex[1], vertex[0]);            if(1.5 * vertex[1].x < tmp) {                if(1.5 * vertex[0].x >= tmp) {                    r = tmp, h = GetY(vertex[1], vertex[0]);                } else {                    r = 1.5 * vertex[0].x, h = 3.0 * vertex[0].y;                }            } else {                if(1.5 * vertex[0].x > tmp) {                    r = 1.5 * vertex[1].x, h = 3.0 * vertex[1].y;                } else {                    double r1 = vertex[0].x * 1.5, h1 = vertex[0].y * 3;                    double r2 = vertex[1].x * 1.5, h2 = vertex[1].y * 3;                    if(r1 * r1 * h1 < r2 * r2 * h2) {                        r = r1, h = h1;                    } else {                        r = r2, h = h2;                    }                }            }            printf("%.3lf %.3lf\n", h, r);            continue;        }        vertex[m] = vertex[0];        double MaxX = -1.0, MaxY = -1.0;        int idx, idy;        for(int i = 0; i < m; ++i) {            if(sgn(vertex[i].x - MaxX) > 0) {                MaxX = vertex[i].x;                idx = i;            }            if(sgn(vertex[i].y - MaxY) > 0) {                MaxY = vertex[i].y;                idy = i;            }        }        if(idx == idy) {            printf("%.3lf %.3lf\n", vertex[idx].y * 3.0, vertex[idx].x * 1.5);            continue;        }        double pre, cur, r2, h2, v2;        pre = GetX(vertex[(idx + 1) % m], vertex[idx]);        if(1.5 * vertex[idx].x <= pre) {            r = 1.5 * vertex[idx].x, h = 3.0 * vertex[idx].y;            v = r * r * h;        }        v = 1e20;        for(int i = idx + 1; i < idy; ++i) {            if(sgn(vertex[i].y - vertex[i - 1].y) < 0) continue;            pre = GetX(vertex[i], vertex[i - 1]);            cur = GetX(vertex[i + 1], vertex[i]);            if(1.5 * vertex[i].x < pre) {                r2 = pre, h2 = GetY(vertex[i], Point(pre, 0));            } else if(1.5 * vertex[i].x < cur) {                r2 = 1.5 * vertex[i].x, h2 = 3 * vertex[i].y;            }else {                r2 = cur, h2 = GetY(vertex[i], Point(cur, 0));            }            v2 = r2 * r2 * h2;            if(sgn(v2 - v) < 0) {                v = v2, r = r2, h = h2;            }        }        pre = GetX(vertex[idy], vertex[idy - 1]);        if(1.5 * vertex[idy].x >= pre) {            r2 = 1.5 * vertex[idy].x, h2 = 3.0 * vertex[idy].y;            v2 = r2 * r2 * h2;            if(sgn(v2 - v) < 0) {                v = v2, r = r2, h = h2;            }        } else {            r2 = pre, h2 = GetY(vertex[idy], Point(pre, 0));            v2 = r2 * r2 * h2;            if(sgn(v2 - v) < 0) {                v = v2, r = r2, h = h2;            }        }        printf("%.3lf %.3lf\n", h, r);    }    return 0;}