【Vijos 1998】【SDOI 2016】平凡的骰子

https://vijos.org/p/1998
三維計算幾何。
須要混合積求四面體體積；
四面體剖分後合併帶權重心求總重心；
四面體重心的橫縱座標是四個頂點的橫縱座標的平均數；
三維差積求平面的法向量；
點積求法向量夾角（二面角）
這些知識就能夠了AC此題了。
時間複雜度\(O(nf)\)，注意\(n,f\leq 100\)，題面描述有誤。

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;

const int N = 103;
const double Pi = acos(-1);

struct Point {
    double x, y, z;
    Point(double _x = 0, double _y = 0, double _z = 0) : x(_x), y(_y), z(_z) {}
    Point operator / (const double &A) const {
        return Point(x / A, y / A, z / A);
    }
    Point operator + (const Point &A) const {
        return Point(x + A.x, y + A.y, z + A.z);
    }
    Point operator - (const Point &A) const {
        return Point(x - A.x, y - A.y, z - A.z);
    }
    double operator * (const Point &A) const {
        return x * A.x + y * A.y + z * A.z;
    }
    Point operator ^ (const Point &A) const {
        return Point(y * A.z - z * A.y, z * A.x - x * A.z, x * A.y - y * A.x);
    }
    double len() {
        return sqrt(x * x + y * y + z * z);
    }
} P[N], H[N * N];

double val[N * N];
int n, F[N][N], m, Htot = 0;

int main() {
    scanf("%d%d", &n, &m);
    for (int i = 1; i <= n; ++i) scanf("%lf%lf%lf", &P[i].x, &P[i].y, &P[i].z);
    for (int i = 1; i <= m; ++i) {
        scanf("%d", &F[i][0]);
        for (int j = 1; j <= F[i][0]; ++j)
            scanf("%d", &F[i][j]);
    }
    
    Point u = P[1];
    for (int i = 1; i <= m; ++i) {
        Point u2 = P[F[i][1]], v1, v2;
        for (int j = 2; j < F[i][0]; ++j) {
            v1 = P[F[i][j]]; v2 = P[F[i][j + 1]];
            H[++Htot] = (u + u2 + v1 + v2) / 4;
            val[Htot] = fabs(((v1 - u2) ^ (v2 - u2)) * (u - u2));
        }
    }
    
    double valtot = 0;
    u = Point(0, 0, 0);
    
    for (int i = 1; i <= Htot; ++i) {
        valtot += val[i];
        u = u + Point(H[i].x * val[i], H[i].y * val[i], H[i].z * val[i]);
    }
    
    u = u / valtot;
    
    for (int i = 1; i <= m; ++i) {
        double ans = 0, co;
        Point u1, u2, u3;
        for (int j = 1, s1, s2; j <= F[i][0]; ++j) {
            s1 = j + 1; if (s1 > F[i][0]) s1 = 1;
            s2 = s1 + 1; if (s2 > F[i][0]) s2 = 1;
            u1 = P[F[i][j]] - u;
            u2 = P[F[i][s1]] - u;
            u3 = P[F[i][s2]] - u;
            u1 = (u1 ^ u2);
            u3 = (u3 ^ u2);
            co = u1 * u3 / u1.len() / u3.len();
            ans += acos(co);
        }
        ans -= (F[i][0] - 2) * Pi;
        printf("%.7lf\n", ans / Pi / 4);
    }
    
    return 0;
}