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

最新推荐文章于 2020-09-23 12:51:45 发布

ramay7

最新推荐文章于 2020-09-23 12:51:45 发布

阅读量2k

点赞数

分类专栏： UVa Online Judge +++三分、模拟退火文章标签： UVALive 三分凸性函数

本文链接：https://blog.csdn.net/Ramay7/article/details/51975712

版权

UVa Online Judge 同时被 2 个专栏收录

38 篇文章 0 订阅

订阅专栏

+++三分、模拟退火

2 篇文章 0 订阅

订阅专栏

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

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

(x, y, z) \to (x 2 + y 2 - - - - - - \sqrt, z)

$(x,y,z)\rightarrow (\sqrt{x^2+y^2},z)$
这样子一来我们就把所有的点限制在第一象限内。
我们来考虑当第一象限内的一个点

(x,y) $(x,y)$ 恰好在圆锥面上时的各个变量之间的关系，设此时圆锥的底面圆半径和高分别为

r和h $r和h$ 。通过相似三角形可以得到：

h - z x 2 + y 2 - - - - - - \sqrt = h r

$\frac{h-z}{\sqrt{x^2+y^2}}=\frac{h}{r}$
化简一下可得:

h = z * r r - x 2 + y 2 - - - - - - \sqrt

$h=\frac{z*r}{r-\sqrt{x^2+y^2}}$
我们带进圆锥的体积公式可得：

v = 1 3 * π * h * r 2 = 1 3 * π * z * r 3 r - x 2 + y 2 - - - - - - \sqrt

$v =\frac{1}{3}*\pi*h*r^2\\ =\frac{1}{3}*\pi*z*\frac{r^3}{r-\sqrt{x^2+y^2}}$
我们对

v $v$ 求两次导可得（抛开系数

13∗π∗z $\frac{1}{3}*\pi*z$ 不看，并令

t=x2+y2−−−−−−√ $t=\sqrt{x^2+y^2}$ ）：

v' = r 2 * ( 2 * r - 3 * t ) ( r - t ) 2

$v'=\frac{r^{2}*(2*r-3*t)}{{(r-t)}^2}$

v'' = 2 * r * ( r 2 + 3 * t 2 - 3 * r * t ) ( r - t ) 3 = 2 r * [ ( r - 3 2 t ) 2 + 1 4 t 2 ] ( r - t ) 3

$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′′ $v''$ 其实是一个恒大于

0 $0$ 的函数，所以函数

v′ $v'$ 是单调递增的。同时我们可以发现函数

v′ $v'$ 是先负后正的。那么函数

v $v$ 关于

r $r$ 就是先增后减的，是个凸函数。

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

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

h = z * r r - x 2 + y 2 - - - - - - \sqrt

$h=\frac{z*r}{r-\sqrt{x^2+y^2}}$
对于每个点如果

r $r$ 确定了，那么最小体积圆锥的

h $h$ 也就确定了，为了满足所有点都覆盖，需要取

hmax $h_{max}$ 。比较

r=mid $r=mid$ 和

r=midmid $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 <script type="math/tex" id="MathJax-Element-15793">l</script>和边与横轴交点，对每个顶点求能覆盖的最小圆锥，也可已解决的。
但是写起来太麻烦了。。。。而且边界情况比较多。
下面的两个代码都是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;
}

ramay7

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
UVALive 4986 Dome of Circus(三分、凸包、凸性函数)

题目链接； UVALive 4986 Dome of Circus 题意：在空间中给nn个点的坐标，并且每个点都是在zz轴的正半平面内，以原点(0,0,0)(0,0,0)为底面圆圆心，画一个圆锥，要求圆锥要包含所有的点，点可以在圆锥面上，求最小的圆锥体积的底面圆半径和高。数据范围：1≤n≤10000,并且所有点的坐标的绝对值都不超过100001 ≤ n ≤ 10000,并且所有点的坐
复制链接

扫一扫