【bzoj4603】【SDOI2016】【平凡的骰子】【计算几何】

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的面积为A
rea(ABC)=alpha+beta+gamma-pi，其中alpha，beta和gamma分别对应了三个二面角的大小。如下图所示。
 
我们保证：每一面着地的时候，重心的垂心都恰好在这一面内。也就是说不会出现摆不稳的情况。

Input

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

Output

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

Sample Input

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

Sample Output

0.1666667

0.1666667

0.1666667

0.1666667

0.1666667

0.1666667

题解：

可以发现题目分为两步.

1.求重心.2.算球面多边形面积.

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

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

求四面体的体积可以用一个点引出的三个向量的混合积乘1/6.

然后重心就求完了.

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

求二面角可以用法向量.

代码还是不难写的。。。

代码：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#define pi 4.0*atan(1.0)
#define N 110
using namespace std;
struct Node{double x,y,z;}node
,cen;
struct Surface{
int num;
Node node
;
}surf
;
int n,m,x;
Node operator+(Node a,Node b){
Node c;
c.x=a.x+b.x;
c.y=a.y+b.y;
c.z=a.z+b.z;
return c;
}
Node operator-(Node a,Node b){
Node c;
c.x=a.x-b.x;
c.y=a.y-b.y;
c.z=a.z-b.z;
return c;
}
Node operator*(Node a,double b){
Node c;
c.x=a.x*b;
c.y=a.y*b;
c.z=a.z*b;
return c;
}
Node operator/(Node a,double b){
Node c;
c.x=a.x/b;
c.y=a.y/b;
c.z=a.z/b;
return c;
}
double operator&(Node a,Node b){
double ans=0.0;
ans+=a.x*b.x;
ans+=a.y*b.y;
ans+=a.z*b.z;
return ans;
}
Node operator|(Node a,Node b){
Node c;
c.x=a.y*b.z-a.z*b.y;
c.y=a.z*b.x-a.x*b.z;
c.z=a.x*b.y-a.y*b.x;
return c;
}
double dis(Node a){
double ans;
ans=a.x*a.x+a.y*a.y+a.z*a.z;
ans=sqrt(ans);
return ans;
}
void init(){
scanf("%d%d",&n,&m);
for (int i=1;i<=n;i++) scanf("%lf%lf%lf",&node[i].x,&node[i].y,&node[i].z);
for (int i=1;i<=m;i++){
scanf("%d",&surf[i].num);
for (int j=0;j<surf[i].num;j++){
scanf("%d",&x);
surf[i].node[j]=node[x];
}
}
}
void getcenter(){
for (int i=1;i<=n;i++)
cen=cen+node[i];
cen=cen/(double)n;
}
double getv(Node p,Node x,Node y,Node z){
x=x-p;y=y-p;z=z-p;
double ans=(x|y)&z;
ans=ans/6.0;
if (ans<0) ans=-ans;
return ans;
}
Node getc(Node p,Node x,Node y,Node z){
Node c;
c=p+x+y+z;
c=c/4.0;
return c;
}
void changecenter(){
double sv=0.0;Node t;
t.x=t.y=t.z=0.0;
for (int i=1;i<=m;i++){
int cc=surf[i].num;
for (int j=1;j<=cc-2;j++){
double v=getv(cen,surf[i].node[0],surf[i].node[j],surf[i].node[j+1]);
Node tt=getc(cen,surf[i].node[0],surf[i].node[j],surf[i].node[j+1]);
tt=tt*v;sv+=v;t=t+tt;
}
}
cen=t/sv;
}
double getangle(Node p,Node x,Node y,Node z){
x=x-p;y=y-p;z=z-p;
x=x/dis(x);y=y/dis(x);z=z/dis(x);
Node a=x|y,b=x|z;
double ans=(a&b)/dis(a)/dis(b);
ans=acos(ans);
return ans;
}
double getsq(Surface &x,Node y){
double ans=-(x.num-2)*pi;
for (int i=0;i<x.num;i++)
ans+=getangle(y,x.node[i],x.node[(i+1)%x.num],x.node[(i-1+x.num)%x.num]);
return ans;
}
void getans(){
for (int i=1;i<=m;i++){
double ans=getsq(surf[i],cen);
ans/=4.0*pi;
printf("%.7lf\n",ans);
}
}
int main(){
init();
getcenter();
changecenter();
getans();
}