【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[N],cen;
struct Surface{
  int num;
  Node node[N];
}surf[N];
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();
}