[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
樣例輸出1
0.1666667
0.1666667
0.1666667
0.1666667
0.1666667
0.1666667
限制
首先存在20%的數據，骰子為長方體。
其次存在20%的數據，骰子為四面體。
余下的數據中有30%的數據，每一面都是三角形。
對於100%的數據，4<=n<=100
來源
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\n",ans/4/pi); } return 0; }
[Sdoi2016]平凡的骰子