[LOJ 2070] 「SDOI2016」平凡的骰子
【题目链接】
【题解】
原题求的是球面面积
可以理解为首先求多面体重心，然后算球面多边形的面积
求重心需要将多面体进行四面体剖分，从而计算出每一个四面体的重心和体积，加权平均即为整个多面体的重心
四面体体积可以用一个点引出的三条向量的积乘 \(\frac 1 6\)
四面体重心坐标是四个顶点坐标平均数
根据题目提示，球面多边形面积为三个二面角之和减去 \(\pi\)，那么我们需要求二面角
先求出法向量，然后点积求向量二面角
【代码】
// Copyright lzt
#include<stdio.h>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<vector>
#include<map>
#include<set>
#include<cmath>
#include<iostream>
#include<queue>
#include<string>
#include<ctime>
using namespace std;
typedef long long ll;
typedef pair<int, int> pii;
typedef long double ld;
typedef unsigned long long ull;
typedef pair<long long, long long> pll;
#define fi first
#define se second
#define pb push_back
#define mp make_pair
#define rep(i, j, k) for (register int i = (int)(j); i <= (int)(k); i++)
#define rrep(i, j, k) for (register int i = (int)(j); i >= (int)(k); i--)
#define Debug(...) fprintf(stderr, __VA_ARGS__)
inline ll read() {
ll x = 0, f = 1;
char ch = getchar();
while (ch < '0' || ch > '9') {
if (ch == '-') f = -1;
ch = getchar();
}
while (ch <= '9' && ch >= '0') {
x = 10 * x + ch - '0';
ch = getchar();
}
return x * f;
}
#define enter putchar('\n')
#define space putchar(' ')
#define MAXN 1000005
#define mo 999999137
typedef long long int64;
typedef double db;
template <class T>
void read(T &res) {
res = 0;
T f = 1;
char c = getchar();
while (c < '0' || c > '9') {
if (c == '-')
f = -1;
c = getchar();
}
while (c >= '0' && c <= '9') {
res = res * 10 + c - '0';
c = getchar();
}
res *= f;
}
template <class T>
void out(T x) {
if (x < 0) {
x = -x;
putchar('-');
}
if (x >= 10)
out(x / 10);
putchar('0' + x % 10);
}
const db PI = acos(-1.0);
struct Point {
db x, y, z;
Point() {}
Point(db _x, db _y, db _z) {
x = _x;
y = _y;
z = _z;
}
friend Point operator+(const Point &a, const Point &b) { return Point(a.x + b.x, a.y + b.y, a.z + b.z); }
friend Point operator-(const Point &a, const Point &b) { return Point(a.x - b.x, a.y - b.y, a.z - b.z); }
friend Point operator*(const Point &a, const db &d) { return Point(a.x * d, a.y * d, a.z * d); }
friend Point operator/(const Point &a, const db &d) { return Point(a.x / d, a.y / d, a.z / d); }
friend Point operator*(const Point &a, const Point &b) {
return Point(a.y * b.z - a.z * b.y, -a.x * b.z + a.z * b.x, a.x * b.y - a.y * b.x);
}
friend db dot(const Point &a, const Point &b) { return a.x * b.x + a.y * b.y + a.z * b.z; }
Point operator-=(const Point &b) { return *this = *this - b; }
Point operator+=(const Point &b) { return *this = *this + b; }
Point operator/=(const db &d) { return *this = *this / d; }
Point operator*=(const db &d) { return *this = *this * d; }
db norm() { return sqrt(x * x + y * y + z * z); }
} P[55], G;
vector<Point> S[85];
int N, F;
Point GetG(Point p, Point a, Point b, Point c) { return (p + a + b + c) / 4.0; }
db GetV(Point p, Point a, Point b, Point c) {
a -= p;
b -= p;
c -= p;
db res = abs(dot(a, b * c));
res /= 6.0;
return res;
}
Point CalcG() {
Point t = Point(0.0, 0.0, 0.0);
db sv = 0.0;
for (int i = 1; i <= F; ++i) {
int s = S[i].size();
for (int j = 0; j <= s - 1; ++j) {
Point tmp = GetG(P[1], S[i][j], S[i][(j + 1) % s], S[i][(j + 2) % s]);
db v = GetV(P[1], S[i][j], S[i][(j + 1) % s], S[i][(j + 2) % s]);
sv += v;
t += tmp * v;
}
}
t /= sv;
return t;
}
db CalcTangle(Point p, Point x, Point y, Point z) {
x -= p;
y -= p;
z -= p;
return acos(dot(x * y, x * z) / (x * y).norm() / (x * z).norm());
}
void Init() {
read(N);
read(F);
db x, y, z;
for (int i = 1; i <= N; ++i) {
scanf("%lf%lf%lf", &x, &y, &z);
P[i] = Point(x, y, z);
}
int k, a;
for (int i = 1; i <= F; ++i) {
read(k);
for (int j = 1; j <= k; ++j) {
read(a);
S[i].pb(P[a]);
}
}
}
void Solve() {
Point G = CalcG();
for (int i = 1; i <= F; ++i) {
int s = S[i].size();
db x = -(s - 2) * PI;
for (int j = 0; j < s; ++j) {
x += CalcTangle(G, S[i][j], S[i][(j + 1) % s], S[i][(j - 1 + s) % s]);
}
printf("%.7lf\n", x / (4 * PI));
}
}
int main() {
Init();
Solve();
return 0;
}