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

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

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

$$(x,y,z)\rightarrow (\sqrt{x^2+y^2},z)$$
这样子一来我们就把所有的点限制在第一象限内。
我们来考虑当第一象限内的一个点$(x,y)$恰好在圆锥面上时的各个变量之间的关系，设此时圆锥的底面圆半径和高分别为$r和h$。通过相似三角形可以得到：
$$\frac{h-z}{\sqrt{x^2+y^2}}=\frac{h}{r}$$
化简一下可得:
$$h=\frac{z*r}{r-\sqrt{x^2+y^2}}$$
我们带进圆锥的体积公式可得：
$$v =\frac{1}{3}*\pi*h*r^2\=\frac{1}{3}*\pi*z*\frac{r^3}{r-\sqrt{x^2+y^2}}$$
我们对$v$求两次导可得（抛开系数$\frac{1}{3}*\pi*z$不看，并令$t=\sqrt{x^2+y^2}$）：
$$v‘=\frac{r^{2}*(2*r-3*t)}{{(r-t)}^2}$$
$$v‘‘=\frac{2*r*(r^2+3*t^2-3*r*t)}{{(r-t)}^3}\=\frac{2r*[(r-\frac{3}{2}t)^2+\frac{1}{4}t^2]}{(r-t)^3}$$
于是我们可以发现$v‘‘$其实是一个恒大于$0$的函数，所以函数$v‘$是单调递增的。同时我们可以发现函数$v‘$是先负后正的。那么函数$v$关于$r$就是先增后减的，是个凸函数。

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

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

$$h=\frac{z*r}{r-\sqrt{x^2+y^2}}$$
对于每个点如果$r$确定了，那么最小体积圆锥的$h$也就确定了，为了满足所有点都覆盖，需要取$h_{max}$。比较$r=mid$和$r=midmid$时的体积，取最小即可。

但是我们也可以换个思路。回到最初的问题，我们是想在第一象限内画一条斜率为负的直线，并且假设和横轴、纵轴分别交于$(r,0),(0,h)$，如何使得$h*r^2$最小？
在看看上面的推导根据$v$的一阶导数可知当$v‘=0$时，可以使得$v$最小，也就是$r=\frac{3}{2} * t,其中t=\sqrt{x^2+y^2}$,那么也就是说对于每个点如果单独考虑的话，我们可以直接取$r=\frac{3}{2} * t,其中t=\sqrt{x^2+y^2}，h=\frac{z*r}{r-\sqrt{x^2+y^2}}$使得这个圆锥体积最小，设这时直线为$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;
}