【凸包】【三分】Gym - 101309D - Dome of Circus

容易发现，圆锥体积和点的具体x、y坐标无关，只与其到z轴的距离sqrt(x*x+y*y)有关。

于是将这些三维的点都投射到二维的xOy平面的第二象限(sqrt(x*x+y*y),z)，求个上凸壳，然后在每一点处，圆锥的母线的斜率的取值范围就确定了，发现这个圆锥的体积关于圆锥母线的函数是单峰的，可以三分。

于是枚举凸壳上每一个点，做个三分就行了。

#include<cstdio>
#include<cmath>
#include<algorithm>
#define EPS 0.00000001
using namespace std;
struct Point{
double x,y;
};
typedef Point Vector;
Vector operator - (const Point &a,const Point &b){
return (Vector){a.x-b.x,a.y-b.y};
}
double Cross(const Vector &a,const Vector &b){
return a.x*b.y-a.y*b.x;
}
bool cmp(const Point &a,const Point &b){
return fabs(a.x-b.x)>=EPS ? a.x<b.x : a.y<b.y;
}
int n,e;
Point ps[10010],qs[10010];
double V=10000000000000.0,R,H;
double sqr(double x){
return x*x;
}
double f(int K,double x){
return sqr(qs[K].y/x-qs[K].x)*(-qs[K].x*x+qs[K].y);
}
int main()
{
freopen("dome.in","r",stdin);
freopen("dome.out","w",stdout);
double X,Y,Z,maxZ=0,maxXY=0;
scanf("%d",&n);
for(int i=0;i<n;++i){
scanf("%lf%lf%lf",&X,&Y,&Z);
ps[i]=(Point){-sqrt(X*X+Y*Y),Z};
maxZ=max(maxZ,Z);
maxXY=max(maxXY,sqrt(X*X+Y*Y));
}
ps[n++]=(Point){0,maxZ};
ps[n++]=(Point){-maxXY,0};
sort(ps,ps+n,cmp);
for(int i=n-1;i>=0;--i){
while(e>1 && Cross(qs[e-1]-qs[e-2],ps[i]-qs[e-1])<EPS){
--e;
}
qs[e++]=ps[i];
}
for(int i=1;i<e-1;++i){
double l=(qs[i-1].y-qs[i].y)/(qs[i-1].x-qs[i].x),r;
if(fabs(l)<EPS){
l+=EPS;
}
if(fabs(qs[i].x-qs[i+1].x)>=EPS){
r=(qs[i].y-qs[i+1].y)/(qs[i].x-qs[i+1].x);
}
else{
r=10000000000000.0;
}
while(r-l>EPS){
double m1=(l+(r-l)/3.0);
double m2=(r-(r-l)/3.0);
//			printf("%lf %lf\n",f(i,m1),f(i,m2));
//			puts("");
if(f(i,m1)>f(i,m2)){
l=m1;
}
else{
r=m2;
}
//			printf("%lf %lf\n",l,r);
}
double fl=f(i,l);
if(fl<V){
V=fl;
R=qs[i].y/l-qs[i].x;
H=-qs[i].x*l+qs[i].y;
}
}
printf("%.3f %.3f\n",H,R);
return 0;
}