bzoj4603: [Sdoi2016]平凡的骰子 | jeffswt: the ac moments

Description

这是一枚平凡的骰子. 它是一个均质凸多面体, 表面有 \(n\) 个端点, 有 \(f\) 个面, 每一面是一个凸多边形, 且任意两面不共面. 将这枚骰子抛向空中, 骰子落地的时候不会发生二次弹跳 (这是一种非常理想的情况). 你希望知道最终每一面着地的概率. 每一面着地的概率可以用如下的方法计算: 我们假设 \(O\) 为骰子的重心, 并以 \(O\) 为球心, 做半径为 \(1\) 的单位球面 (记为 \(S\)) . 我们知道 \(S\) 的表面积即单位球的表面积, 为 \(4\pi\), 这里 \(\pi\) 为圆周率. 对于骰子的某一面 \(C\) 来说, 球面 \(S\) 上存在一块区域 \(T\) 满足: 当下落时若骰子所受重力方向与 \(S\) 的交点落在 \(T\) 中, 则 \(C\) 就是最终着地的一面. 那么 \(C\) 着地的概率为区域 \(T\) 的面积除以 \(4\pi\). 为了能更好地辅助计算球面上一块区域的面积, 我们给出单位球面 \(S\) 上三角形的面积计算公式. 考虑单位球面 \(S\) 上的三个两两相交的大圆, 交点依次为 \(A, B, C\). 则曲面三角形 \(ABC\) 的面积为 \(S_{ABC} = \alpha+\beta+\gamma-\pi\), 其中 \(\alpha, \beta, \gamma\) 分别对应了三个二面角的大小.

我们保证: 每一面着地的时候, 重心的垂心都恰好在这一面内. 也就是说不会出现摆不稳的情况.

Input

第一行输入两个整数, 分别表示端点总数 \(n\) 与表面总数 \(f\), 分别从 \(1\) 开始编号. 之后 \(n\) 行, 每行有三个浮点数 \(x, y, z\), 给出了每一个端点的坐标. 之后 \(f\) 行依次描述了每一块表面, 首先给出不小于 \(3\) 的整数 \(d\), 表示这一面的端点个数, 之后 \(d\) 个整数按照逆时针方向 (视角在骰子的外面) 给出了每一个端点的编号.

Output

输出 \(f\) 行, 第 \(i\) 行有一个浮点数, 表示第 \(i\) 个面着地的概率. 本题中您的输出应该保留距离答案最近的 \(7\) 位小数, 即在需要保留 \(7\) 位小数的前提之下与标准答案最接近. 数据保证可以避免对小数点后第八位四舍五入后产生的精度误差.

Sample Input

Sample Output

Data Range

对于所有数据, 满足:\(4 \leq n \leq 50, 4 \leq m \leq 50, 0 \leq |x|, |y| \leq 10000\)

Explanation

参照本神犇题解, 原文链接:http://blog.csdn.net/sunshinezff/article/details/51463895

可以发现题目分为两步:

求重心

算球面多边形面积

求重心我们可以对整个多面体进行四面体剖分.

求出每个四面体的重心, 然后再按它们的体积加权平均.

求四面体的体积可以用一个点引出的三个向量的混合积乘 \(\frac{1}{6}\),

然后重心就求完了.

算面积直接按题目说的方法算就可以.

求二面角可以用法向量.

代码还是不难写的......

Source Code


#include <iostream>
#include <cstdlib>
#include <cstdio>
#include <cmath>

using namespace std;
typedef long long lli;
typedef double llf;
const int maxn = 110;
const llf pi = 3.1415926535897932384626;

struct node {
    llf x, y, z;
    node& operator += (const node& b) {
        x += b.x, y += b.y, z += b.z; return *this; }
    node& operator -= (const node& b) {
        x -= b.x, y -= b.y, z -= b.z; return *this; }
    node& operator *= (const llf& b) {
        x *= b, y *= b, z *= b; return *this; }
    node& operator /= (const llf& b) {
        x /= b, y /= b, z /= b; return *this; }
    node(llf _1, llf _2, llf _3) {
        x = _1, y = _2, z = _3; return ; }
    node(void) {
        x = y = z = 0.0; return ; }
    llf dist(void) {
        return sqrt(x*x + y*y + z*z); }
};
node operator + (const node& a, const node& b) {
    node c = a; c += b; return c; }
node operator - (const node& a, const node& b) {
    node c = a; c -= b; return c; }
node operator * (const node& a, const llf& b) {
    node c = a; c *= b; return c; }
node operator / (const node& a, const llf& b) {
    node c = a; c /= b; return c; }
llf operator * (const node& a, const node& b) {
    return a.x*b.x + a.y*b.y + a.z*b.z; }
node operator ^ (const node& a, const node& b) {
    return node(a.y*b.z-a.z*b.y, a.z*b.x-a.x*b.z, a.x*b.y-a.y*b.x); }
node get_centre(const node& a, const node& b, const node& c, const node& d) {
    return (a + b + c + d) / 4.0; }
llf get_volume(node p, node a, node b, node c) {
    a -= p, b -= p, c -= p;
    return fabs(((a^b) * c) / 6.0); }
llf get_angle_rad(node p, node a, node b, node c) {
    a -= p, b -= p, c -= p;
    a /= a.dist(), b /= b.dist(), c /= c.dist();
    node x = a^b, y = a^c;
    return acos(x*y / x.dist() / y.dist()); }

struct surface {
    int n;
    node nd[maxn];
};
llf get_curve_area(const surface& sf, const node& cent)
{
    llf res = -(sf.n - 2) * pi;
    int n = sf.n;
    for (int i = 0; i < n; i++)
        res += get_angle_rad(cent, sf.nd[i], sf.nd[(i+1)%n], sf.nd[(i-1+n)%n]);
    return res;
}

int n, m;
node nd[maxn], cent;
surface sf[maxn];

int main(int argc, char** argv)
{
    scanf("%d%d", &n, &m);
    for (int i = 1; i <= n; i++)
        scanf("%lf%lf%lf", &nd[i].x, &nd[i].y, &nd[i].z);
    for (int i = 1; i <= m; i++) {
        scanf("%d", &sf[i].n);
        for (int j = 0, x; j < sf[i].n; j++) {
            scanf("%d", &x);
            sf[i].nd[j] = nd[x];
        }
    }
    // Calculating the centre of the object
    for (int i = 1; i <= n; i++)
        cent += nd[i];
    cent /= n;
    // Adjusting the centre of the object from pre-calculations
    node f_cent;
    llf f_vol;
    for (int i = 1; i <= m; i++) {
        for (int j = 1; j < sf[i].n - 1; j++) {
            llf vol = get_volume(cent, sf[i].nd[0], sf[i].nd[j], sf[i].nd[j+1]);
            node t = get_centre(cent, sf[i].nd[0], sf[i].nd[j], sf[i].nd[j+1]);
            t *= vol; f_vol += vol; f_cent += t;
        }
    }
    cent = f_cent / f_vol;
    // Calculating result
    for (int i = 1; i <= m; i++) {
        llf res = get_curve_area(sf[i], cent);
        res /= 4.0 * pi;
        printf("%.7lf\n", res);
    }
    return 0;
}