题目大意
给定一个凸多面体,求其抛起后每个面落地的概率。
题解
这道题其实主要难点在于找重心。
凸多边形重心的计算方案应该就是三角剖分,每个三角形的重心就是三个点坐标的平均值。然后把这些重心按照所在三角形的面积加权平均。
类比到凸多面体,我们对其进行四面体剖分,每个四面体的重心也就是四个顶点坐标的平均值,重心加权平均出来即可。四面体体积可以用行列式/向量混合积求出。
至于为什么,我也不知道
接下来就是求二面角了,这个玩意儿可以直接求两个平面法线夹角。于是就做完了。
#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;
}