【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;
}