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;
}