HDU Dome of Circus

#include <iostream>
#include <cstdio>
#include <queue>
#include <cstring>
#include <cmath>
#include <vector>
#include <set>
#include <stack>
#include <algorithm>
//#include "myAlgorithm.h"

#define MAX 10005
#define INF (1e8 + 5)
#define eps 1e-7
#define Rep(s, e) for( int i = s; i <= e; i++)
#define Cep(e, s) for( int i = e; i >= s; i --)
#define PI acos(-1.0)
using namespace std;
int T;
int n, tn;
double h, r, minVol;
bool isFirst;
struct Point{
double x, y;
void set(double x, double y){this->x = x; this->y = y;}
Point& operator + (const Point & p){
Point t;
t.set(p.x + x, p.y + y);
return t;
}
Point& operator - (const Point & p){
Point t;
t.set( x - p.x, y - p.y);
return t;
}
bool operator < (const Point &p)const{
if(x != p.x)
return x > p.x;
return y > p.y;
}
void show(){
cout<<x<<" "<<y<<endl;
}
}p[MAX], tubao[MAX], S, E;
double croV(const Point &vs,const Point &ve)
{
return (vs.x * ve.y - vs.y * ve.x);
}
void gtubao(int node)
{
tubao[tn++] = p[node];
if(node == n - 1)return ;
Point vs, ve, pt;

bool isFind = 0;
int newNode = node;
for(int i = node + 1; i < n; i++){
if(p[i].x != p[newNode].x && p[i].y > p[newNode].y ){
if(!isFind){
pt = p[i];
isFind = 1;
newNode = i;
}else {
vs = p[i] - p[node];
ve = p[i] - pt;
double crp = croV(vs, ve);
if(crp >= 0 ){
pt = p[i];
newNode = i;
}
}
}
}

if(newNode != node){
gtubao(newNode);
}

}
void init()
{
tn = 0;
cin>>n;
double x, y, z;
Rep(0, n - 1){
cin>>x>>y>>z;
x = sqrt(x * x + y * y);
y = z;
p[i].set(x, y);
}sort(p , p + n);
}
double gvol(int node, double angle, double &rt, double &ht){
double l1 = tubao[node].x, hh = tubao[node].y,  l2;
l2 = hh / tan(angle);
rt = l2 + l1;
ht = rt * tan(angle);
return (rt * rt * ht);
}
void makeMin(int node, double sangle, double eangle)///angle not in degree
{
double mid, mmid, vols, vole, rt, ht;
while(eangle -  sangle > eps){
mid = (eangle + sangle)/2; mmid = (mid + eangle)/2;
vols = gvol(node, mid, rt, ht);
vole = gvol(node, mmid, rt, ht);
if(vols < vole){
eangle = mmid;
}else {
sangle = mid;
}
}
if(isFirst || minVol > vole){
minVol = vols;
r = rt;
h = ht;
isFirst = 0;
}
}

void solve(){
cin>>T;
while(T--){
init();
gtubao(0);
isFirst = 1;
//        cout<<tn<<endl;
////        Rep(0, n - 1){
////            p[i].show();
////        }cout<<endl;
////        Rep(0, tn - 1){
////            tubao[i].show();
////        }
///三分
double sangle = 0, eangle;
Cep(tn - 1, 1){
eangle = atan((tubao[i].y - tubao[i - 1].x) / (tubao[i].x - tubao[i - 1].x));
makeMin(i, sangle, eangle);
sangle = eangle;
}
eangle = PI/2;
makeMin(0, sangle, eangle);
printf("%.3f %.3f\n", h, r);
}
}

///-----------------------------抓取一个土包, 然后对凸包顶点进行三分取最小值~~~ ----凸包都O(n^2)了------哭~TLE了-用C++WA--------------

///----------实在debug不出来, 换个写法----------------凸包不用找, 枚举中找出来就行, --------------------
/**
三分用于找凸曲线的极值,
把点投影在x - y上, 线L必然在点集的离散凸包曲线上交于凸包上的某点p(x, y)交x轴于r, y轴于h ,  旋转线L的角度(90, 180),
直觉告诉俄, (h * r * r)是个的值是个曲线,   必然能找到斜率角度degree, 找到min(h * r * r);

三分答案角度：
1. 找到和线L相交的凸包上的点P(x, y), 把线设在过(0, 0)点, 距离线最远的点就是和线相交的点;
2.利用 y = tan(degree) * (x - p.x) + p.y, 计算h和r;

精度太高时间占用高 eps = 1e-10, 太低会wa, 终于畏畏缩缩的AC了 eps = -7~~~
///Accepted 1007 2328 MS 468 KB
*/
int getInsectPoint(double angle){///找和L相交的点
///过p[i]做L‘ : y = tan(angle) * (x - p[i]. x) + p[i].y
///L: y = tan(angle) * x;
double k = tan(angle);
double ans = fabs(p[0].y  - k * p[0].x)/sqrt(1 + k * k);
int num = 0;
Rep(1, n - 1){
double temp = fabs(p[i].y  - k * p[i].x)/sqrt(1 + k * k);
if(temp > ans){
ans = temp;
num = i;
}
}
return num;
}
double gVol(int pos, double angle, double &ht, double &rt){
double k = tan(angle);
ht = -p[pos].x * k + p[pos].y;
rt =  -p[pos].y/k + p[pos].x;
return ht * rt * rt;
}
void makeMin001(){
double l = PI/2, r = PI, lmid, rmid, ht, rt;
int rn, ln;
while(r - l > eps){
lmid = (l + r)/2; rmid = (lmid + r)/2;
ln = getInsectPoint(lmid);
rn = getInsectPoint(rmid);
//cout<<ln<<rn<<endl;
double resl = gVol(ln, lmid, ht, rt), resr = gVol(rn, rmid, ht, rt);
if(resl < resr){
r = rmid;
}else {
l = lmid;
}
}
printf("%.3f %.3f\n", ht, rt);
}
void solve001()
{
cin>>T;
while(T--){
init();
makeMin001();
}
}
int main() {
//freopen("in.tx0t", "w", stdout);
solve001();
return 0;
}
/*
*/