SDOI2016 平凡的骰子（计算几何）

题目链接

题目大意

给定一个凸多面体，求其抛起后每个面落地的概率。

题解

这道题其实主要难点在于找重心。
凸多边形重心的计算方案应该就是三角剖分，每个三角形的重心就是三个点坐标的平均值。然后把这些重心按照所在三角形的面积加权平均。
类比到凸多面体，我们对其进行四面体剖分，每个四面体的重心也就是四个顶点坐标的平均值，重心加权平均出来即可。四面体体积可以用行列式/向量混合积求出。
至于为什么，我也不知道
接下来就是求二面角了，这个玩意儿可以直接求两个平面法线夹角。于是就做完了。

#include <bits/stdc++.h>
using namespace std;
const int MAXR = 10000000;
char _READ_[MAXR], _PRINT_[MAXR];
int _READ_POS_, _PRINT_POS_, _READ_LEN_;
inline char readc() {
#ifndef ONLINE_JUDGE
	return getchar();
#endif
	if (!_READ_POS_) _READ_LEN_ = fread(_READ_, 1, MAXR, stdin);
	char c = _READ_[_READ_POS_++];
	if (_READ_POS_ == MAXR) _READ_POS_ = 0;
	if (_READ_POS_ > _READ_LEN_) return 0;
	return c;
}
template<typename T> inline void read(T &x) {
	x = 0; register int flag = 1, c;
	while (((c = readc()) < '0' || c > '9') && c != '-');
	if (c == '-') flag = -1; else x = c - '0';
	while ((c = readc()) >= '0' && c <= '9') x = x * 10 + c - '0';
	x *= flag;
}
template<typename T1, typename ...T2> inline void read(T1 &a, T2&... x) {
	read(a), read(x...);
}
inline int reads(char *s) {
	register int len = 0, c;
	while (isspace(c = readc()) || !c);
	s[len++] = c;
	while (!isspace(c = readc()) && c) s[len++] = c;
	s[len] = 0;
	return len;
}
inline void ioflush() { fwrite(_PRINT_, 1, _PRINT_POS_, stdout), _PRINT_POS_ = 0; fflush(stdout); }
inline void printc(char c) {
	_PRINT_[_PRINT_POS_++] = c;
	if (_PRINT_POS_ == MAXR) ioflush();
}
inline void prints(char *s) {
	for (int i = 0; s[i]; i++) printc(s[i]);
}
template<typename T> inline void print(T x, char c = '\n') {
	if (x < 0) printc('-'), x = -x;
	if (x) {
		static char sta[20];
		register int tp = 0;
		for (; x; x /= 10) sta[tp++] = x % 10 + '0';
		while (tp > 0) printc(sta[--tp]);
	} else printc('0');
	printc(c);
}
template<typename T1, typename ...T2> inline void print(T1 x, T2... y) {
	print(x, ' '), print(y...);
}
typedef long long ll;

const int MAXN = 105;
const double PI = acos(-1.0);
struct Vec {
	double x, y, z;
	Vec operator+(const Vec &v) const {
		return (Vec) { x + v.x, y + v.y, z + v.z };
	}
	Vec operator-(const Vec &v) const {
		return (Vec) { x - v.x, y - v.y, z - v.z };
	}
	Vec operator*(double d) const {
		return (Vec) { x * d, y * d, z * d };
	}
	Vec operator/(double d) const {
		return (Vec) { x / d, y / d, z / d };
	}
	Vec det(const Vec &v) const {
		return (Vec) { y * v.z - z * v.y, z * v.x - x * v.z, x * v.y - y * v.x };
	}
	double dot(const Vec &v) const {
		return x * v.x + y * v.y + z * v.z;
	}
	double dis() const { return sqrt(this->dot(*this)); }
} po[MAXN];
vector<Vec> pla[MAXN];
int n, m;
int main() {
	scanf("%d%d", &n, &m);
	for (int i = 1; i <= n; i++) scanf("%lf%lf%lf", &po[i].x, &po[i].y, &po[i].z);
	double tot = 0; Vec g = (Vec) { 0, 0, 0 };
	for (int i = 1; i <= m; i++) {
		int t, k; read(k);
		while (k--) { read(t); pla[i].push_back(po[t]); }
		for (int j = 1; j + 1 < (int)pla[i].size(); j++) {
			Vec p = (po[1] + pla[i][0] + pla[i][j] + pla[i][j + 1]) / 4;
			double v = fabs((pla[i][0] - po[1]).det(pla[i][j] - po[1]).dot(pla[i][j + 1] - po[1])) / 6;
			g = g + p * v, tot += v;
		}
	}
	g = g / tot;
	for (int i = 1; i <= m; i++) {
		int t = pla[i].size();
		double area = 0;
		for (int j = 0; j < t; j++) {
			Vec d = pla[i][j] - g;
			Vec a = d.det(pla[i][(j - 1 + t) % t] - g);
			Vec b = d.det(pla[i][(j + 1) % t] - g);
			area += acos(a.dot(b) / (a.dis() * b.dis()));
		}
		printf("%.7lf\n", area / (4 * PI) - (t - 2) * 0.25);
	}
	return 0;
}