題解

用了一堆迷之復雜的結論結果迷之好寫的計算幾何？？？？

好吧，要寫立體幾何了

如果有名詞不懂自己搜吧

首先我們求重心，我們可以求帶權重心，也就是x坐標的話是所有分割的小四面體的x坐標 * 四面體體積的和除以骰子的體積，y，z坐標同理

然後我們把這個骰子四面體剖分，剖分的話就是隨便選在骰子內的一個點，對於骰子的每個面找相鄰的三個點和這個點作為頂點組成的四面體

四面體的重心是四個點三維坐標和除以4

四面體的體積是三維混合積的絕對值除以6

然後對於每個面，我們把它剖分成三角形，發現它們二面角的和就是左右相鄰的兩條邊和重心組成的面二面角的和

求出兩個面二面角的法向量（就是垂直於平面的直線，可以用三維叉積求出來），然後求兩個法向量的夾角，可以求出余弦值然後用反函數

代碼

#include <bits/stdc++.h>
#define enter putchar(‘\n‘)
#define space putchar(‘ ‘)
#define pii pair<int,int>
#define fi first
#define se second
#define mp make_pair
#define MAXN 1000005
#define mo 999999137
#define pb push_back
//#define ivorysi
using namespace std;
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() {
#ifdef ivorysi
    freopen("f1.in","r",stdin);
#endif
    Init();
    Solve();
    return 0;
}
 

【LOJ】#2070. 「SDOI2016」平凡的骰子