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