HDU 4573 Throw the Stones（动态三维凸包）（2013 ACM-ICPC长沙赛区全国邀请赛）

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=4573

Problem Description

　　Remember our childhood? A few naked children throw stones standing on the same position, the one throws farther win the game. Aha, of course, there are some naughty boys who care more about whether they can urinate father.
　　You believe it or not, anyway, I believed. Nowadays, some of the children are smarter than we were, while others may be more naughty.
　　A week ago, I saw several children throw stones. In fact, they are more clever than we were, since the game they played, apparently, is more complex than we did. Maybe you have different points of view, however, you’d better learn about the rules of the game before expressing your views. A group of children take turns to throw stones standing on the same position. After some child throw a stone, the children will draw a convex polyhedron with smallest volume together to enclose all the stones thrown by them. You may assume that the stone is so small as to be abstracted as a point in three-dimensional space. Naively, the children regard the space enclosed by the convex polyhedron as territory under their control. After a child throw his stone, the score he obtains equals the incremental of the volume of their territory.   Unfortunately, the first three child’s score will always be zero. At last, the child with the highest score will win the game, and known as the "King".
　　I think you have accepted my opinion already, for the rules of their throwing stones game are really complicated. But, you also don’t need to be frustrated for it. Now, in order to show you are smarter, maybe you can write a program to help the children point out their "King".

 

Input

　　Input consists of a number of cases. The data of each case appears on a number of input lines, the first of which contains an integer N. The following N lines contain three number (x_i, y_i, z_i) indicating coordinates of the stone thrown by the i-th child. 
Note: 1 <= N <= 10^4, 1 <= i <= N, -10^4 <= x_i , y_i , z_i <= 10^4.

 

Output

　　For each test case, you should output two lines. The first line is "Case #K:", K means the number of the test case. The second line is "i v", i means index of the "King" and v means the score of the "King". If there are more than one "King", output the one throws stone earlier than others.
　　Please round the result to 2 digits after decimal point if necessary.

 

题目大意：给一个三维空间点集依次加入N个点，形成一个三维凸包。问加入哪一个点的时候，凸包的体积增量最大，即前 i 个点组成的凸包体积减去前 i - 1个点组成的凸包体积最大。输出这个点，并输出这个增量。

思路：http://blog.csdn.net/catalyst1314/article/details/9017673

PS：我的代码框架几乎都是照着上面写的。但是怎么看内存都会超出限制，虽然不知道为什么实际上没有（我算错了？），虽说现场没有这种限制。然后不知为何说时间复杂度是O(n*sqrt(n))，我怎么觉得最坏情况下还是O(n^2)的……

 

代码（531MS）：

#include <cstdio>
#include <algorithm>
#include <iostream>
#include <cstring>
#include <cmath>
using namespace std;
typedef long long LL;

const int MAXN = 10001;
const double EPS = 1e-9;

inline int sgn(double x) {
    return (x > EPS) - (x < -EPS);
}

struct Point {
    double x, y, z;
    Point() {}
    Point(double x, double y, double z): x(x), y(y), z(z) {}
    void read() {
        scanf("%lf%lf%lf", &x, &y, &z);
    }
    Point operator - (const Point &rhs) const {
        return Point(x - rhs.x, y - rhs.y, z - rhs.z);
    }
    double operator * (const Point &rhs) const {
        return x * rhs.x + y * rhs.y + z * rhs.z;
    }
};
double length(const Point &a) {
    return sqrt(a * a);
}
Point cross(const Point &a, const Point &b) {
    Point res;
    res.x = a.y * b.z - a.z * b.y;
    res.y = a.z * b.x - a.x * b.z;
    res.z = a.x * b.y - a.y * b.x;
    return res;
}
Point cross(const Point &o, const Point &a, const Point &b) {
    return cross(a - o, b - o);
}
double volume(const Point &a, const Point &b, const Point &c, const Point &d) {
    return cross(c - a , b - a) * (d - a) / 6;
}
Point p[MAXN];

struct Face {
    int a, b, c;
    bool flag;
    Face() {}
    Face(int a, int b, int c, bool flag): a(a), b(b), c(c), flag(flag) {}
    bool can_see(const Point &q) {
        return sgn(volume(p[a], p[b], p[c], q)) > 0;
    }
};
Face fac[MAXN * 10];

struct Convex {
    double diff_vol;
    int cnt, mat[MAXN][MAXN];

    void init() {
        cnt = 0;
        for(int i = 0; i < 4; ++i) {
            Face newf = Face((i + 1) % 4, (i + 2) % 4, (i + 3) % 4, true);
            if(newf.can_see(p[i])) swap(newf.a, newf.c);
            mat[newf.a][newf.b] = mat[newf.b][newf.c] = mat[newf.c][newf.a] = cnt;
            fac[cnt++] = newf;
        }
    }

    void restore(int k, int a, int b) {
        int f = mat[a][b];
        if(fac[f].flag) {
            if(fac[f].can_see(p[k])) dfs(k, f);
            else {
                mat[b][a] = mat[a][k] = mat[k][b] = cnt;
                fac[cnt++] = Face(b, a, k, true);
            }
        }
    }

    void dfs(int k, int f) {
        diff_vol += volume(p[fac[f].a], p[fac[f].b], p[fac[f].c], p[k]);
        fac[f].flag = false;
        restore(k, fac[f].b, fac[f].a);
        restore(k, fac[f].c, fac[f].b);
        restore(k, fac[f].a, fac[f].c);
    }

    double update(int k) {
        diff_vol = 0;
        for(int i = 0; i < cnt; ++i) {
            if(!fac[i].flag || !fac[i].can_see(p[k])) continue;
            dfs(k, i);
            break;
        }
        return diff_vol;
    }

    double vol() {
        double res = 0;
        for(int i = 0; i < cnt; ++i) if(fac[i].flag)
            res -= volume(p[fac[i].a], p[fac[i].b], p[fac[i].c], Point(0, 0, 0));
        return res;
    }
} solver;

int n, kase;

void solve() {
    int king = 0;
    double maxans = 0;
    for(int i = 1, tmp = 1; i < n; ++i) {
        if(tmp == 1) {
            tmp += sgn(length(p[0] - p[i]));
            if(tmp > 1) swap(p[1], p[i]);
        } else if(tmp == 2) {
            tmp += sgn(length(cross(p[0], p[1], p[i])));
            if(tmp > 2) swap(p[2], p[i]);
        } else if(tmp == 3) {
            tmp += (sgn(volume(p[0], p[1], p[2], p[i])) != 0);
            if(tmp > 3) {
                swap(p[3], p[i]);
                solver.init();
                for(int j = 4; j <= i; ++j) solver.update(j);
                king = i, maxans = solver.vol();
            }
        } else {
            double v = solver.update(i);
            if(sgn(v - maxans) > 0) {
                maxans = v;
                king = i;
            }
        }
    }
    printf("%d %.2f\n", king + 1, maxans);
}

int main() {
    while(scanf("%d", &n) != EOF) {
        for(int i = 0; i < n; ++i) p[i].read();
        printf("Case #%d:\n", ++kase);
        solve();
    }
}