[Sdoi2016]平凡的骰子

描述

　　这是一枚平凡的骰子。它是一个均质凸多面体，表面有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 的面积为 Area(ABC)=alpha+beta+gamma-pi，其中 alpha，beta 和 gamma 分别对应了三个二面角的大小。如下图所示。

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

格式

输入格式

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

输出格式

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

样例1

样例输入1

8 6
1 0 0
1 1 0
1 0 1
1 1 1
0 0 0
0 1 0
0 0 1
0 1 1
4 1 2 4 3
4 2 6 8 4
4 6 5 7 8
4 5 1 3 7
4 3 4 8 7
4 1 5 6 2

Copy

样例输出1

0.1666667
0.1666667
0.1666667
0.1666667
0.1666667
0.1666667

Copy

限制

首先存在20%的数据，骰子为长方体。
其次存在20%的数据，骰子为四面体。
余下的数据中有30%的数据，每一面都是三角形。
对于100%的数据，4<=n<=100且4<=m<=100，所有坐标的绝对值都在10000以内。

来源

SDOI 2016 round2 day2

• 三维计算几何。
• 需要混合积求四面体体积；
• 四面体剖分后合并带权重心求总重心；
• 四面体重心的横纵坐标是四个顶点的横纵坐标的平均数；
• 三维差积求平面的法向量；
• 点积求法向量夹角（二面角）
• 这些知识就可以了AC此题了。
• 时间复杂度O(nf)O(nf)，注意n,f≤100n,f≤100，题面描述有误。

#include<cmath>
#include<cstdio>
using namespace std;
inline void read(int &x){
    register char ch=getchar();x=0;
    while(ch<'0'||ch>'9') ch=getchar();
    while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=getchar();
}
const int N=105;
typedef double real;
const real pi=acos(-1);
struct point{
    real x,y,z;
    point(){}
    point(real _x,real _y,real _z):x(_x),y(_y),z(_z){}
    point operator +(const point &a)const{
        return point(x+a.x,y+a.y,z+a.z);
    }
    point operator -(const point &a)const{
        return point(x-a.x,y-a.y,z-a.z);
    }
    point operator ^(const point &a)const{
        return point(y*a.z-z*a.y,z*a.x-x*a.z,x*a.y-y*a.x);
    }
    real operator *(const point &a)const{
        return x*a.x+y*a.y+z*a.z;
    }
    point operator /(const real &a)const{
        return point(x/a,y/a,z/a);
    }
    const real len(){
        return sqrt(x*x+y*y+z*z);
    }
}P[N],H[N*N];
int n,m,Htot,f[N][N];
real val[N*N];
int main(){
    read(n);read(m);
    for(int i=1;i<=n;i++) scanf("%lf%lf%lf",&P[i].x,&P[i].y,&P[i].z);
    for(int i=1;i<=m;i++){
        read(f[i][0]);
        for(int j=1;j<=f[i][0];j++){
            read(f[i][j]);
        }
    }
    point u=P[1];
    for(int i=1;i<=m;i++){
        point u2=P[f[i][1]],v1,v2;
        for(int j=2;j<f[i][0];j++){
            v1=P[f[i][j]];
            v2=P[f[i][j+1]];
            H[++Htot]=(u+u2+v1+v2)/4;//四面体重心 
            val[Htot]=fabs(((u2-v1)^(u2-v2))*(u-u2));//四面体体积 
        }
    }
    u=point(0,0,0);
    real valtot=0;
    for(int i=1;i<=Htot;i++){
        valtot+=val[i];
        u=u+point(H[i].x*val[i],H[i].y*val[i],H[i].z*val[i]);
    }
    u=u/valtot;//球心坐标 
    for(int i=1;i<=m;i++){
        point u1,u2,u3;real co,ans=0;
        for(int j=1,s1,s2;j<=f[i][0];j++){//该平面拆成三角形，计算所有不同夹角 
            s1= j+1;if(s1>f[i][0]) s1=1;
            s2=s1+1;if(s2>f[i][0]) s2=1;
            u1=P[f[i][j]]-u;
            u2=P[f[i][s1]]-u;
            u3=P[f[i][s2]]-u;
            u1=u1^u2;//面(u1,u2)的法向量 
            u3=u3^u2;//面(u2,u3)的法向量 
            co=u1*u3/u1.len()/u3.len();//二面角夹角
            ans+=acos(co); 
        }
        ans-=pi*(f[i][0]-2);
        printf("%.7lf
",ans/4/pi);
    }
    return 0;    
}