湖南省第八届程序设计竞赛L

题目L. 安全区域

 

太空中有n 类卫星（即类型1、类型2…）共m个。对于满足1<=i<=n的任意整数i，类型i的所有卫星共同保卫着包含它们的最小凸多面体（即凸包。可能会退化，即体积为0）。如果空间中一个点被至少k 类卫星保护，我们说这个点属于安全区域。

 

你的任务是计算所有安全区域的总体积。

 

输入格式

输入第一行为数据组数T (T<=25)。每组数据第一行为三个整数n,k和m(1<=k<=n<=5, 4<=m<=50)。以下m 行每行包含三个实数x,y, z，表示一个类型t 的卫星，坐标为(x,y,z)(1<=t<=n, 0<=x,y,z<=10)。每组输入数据后有一个空行。

 

注意：评测数据（而非样例数据）中卫星坐标都是随机生成的。

 

输出格式

对于每组数据，输出安全区域的总体积，四舍五入保留小数点5位。

 

样例输入                                        样例输出

2

2 1 16

1 0 0 0

1 0 0 2

1 0 2 0

1 0 2 2

1 2 0 0

1 2 0 2

1 2 2 0

1 2 2 2

2 1 1 1

2 1 1 3

2 1 3 1

2 1 3 3

2 3 1 1

2 3 1 3

2 3 3 1

2 3 3 3

 

1 1 4

1 0 0 0

1 0 1 0

1 0 0 1

1 1 0 0

15.00000

0.16667

 代码：

// Rujia Liu#include<cstdio>#include<cmath>#include<algorithm>using namespace std;const double eps = 1e-8;int dcmp(double x) {  if(fabs(x) < eps) return 0; else return x < 0 ? -1 : 1;}struct Point3 {  double x, y, z;  Point3(double x=0, double y=0, double z=0):x(x),y(y),z(z) { }};typedef Point3 Vector3;Vector3 operator + (Vector3 A, Vector3 B) {  return Vector3(A.x+B.x, A.y+B.y, A.z+B.z);}Vector3 operator - (Point3 A, Point3 B) {  return Vector3(A.x-B.x, A.y-B.y, A.z-B.z);}Vector3 operator * (Vector3 A, double p) {  return Vector3(A.x*p, A.y*p, A.z*p);}Vector3 operator / (Vector3 A, double p) {  return Vector3(A.x/p, A.y/p, A.z/p);}bool operator < (const Point3& a, const Point3& b) {  if(dcmp(a.x - b.x) != 0) return a.x < b.x;  if(dcmp(a.y - b.y) != 0) return a.y < b.y;  if(dcmp(a.z - b.z) != 0) return a.z < b.z;  return false;}bool operator == (const Point3& a, const Point3& b) {  return dcmp(a.x-b.x) == 0 && dcmp(a.y-b.y) == 0 && dcmp(a.z-b.z) == 0;}Point3 read_point3() {  Point3 p;  scanf("%lf%lf%lf", &p.x, &p.y, &p.z);  return p;}double Dot(Vector3 A, Vector3 B) { return A.x*B.x + A.y*B.y + A.z*B.z; }double Length(Vector3 A) { return sqrt(Dot(A, A)); }double Angle(Vector3 A, Vector3 B) { return acos(Dot(A, B) / Length(A) / Length(B)); }Vector3 Cross(Vector3 A, Vector3 B) { return Vector3(A.y*B.z - A.z*B.y, A.z*B.x - A.x*B.z, A.x*B.y - A.y*B.x); }double Area2(Point3 A, Point3 B, Point3 C) { return Length(Cross(B-A, C-A)); }double Volume6(Point3 A, Point3 B, Point3 C, Point3 D) { return Dot(D-A, Cross(B-A, C-A)); }#include<cstdlib>#include<cstring>#include<queue>using namespace std;double rand01() { return rand() / (double)RAND_MAX; }double randeps() { return (rand01() - 0.5) * eps; }Point3 add_noise(Point3 p) {  return Point3(p.x + randeps(), p.y + randeps(), p.z + randeps());}const int maxn = 2000 + 10;int vis[maxn][maxn];struct Face{  int v[3];  Vector3 normal(Point3 *P) const {    return Cross(P[v[1]]-P[v[0]], P[v[2]]-P[v[0]]);  }  int cansee(Point3* P, int i) const {    return Dot(P[i]-P[v[0]], normal(P)) > 0 ? 1 : 0;  }  double Area2(Point3* P) const {    return ::Area2(P[v[0]], P[v[1]], P[v[2]]);  }  double Volume6(Point3* P, const Point3& C) const {    return ::Volume6(P[v[0]], P[v[1]], P[v[2]], C);  }  double distance(Point3* P, const Point3& C) const {    return -Volume6(P, C) / Area2(P);  }  Point3 centroid(Point3* P, const Point3& C) const {   return (C+P[v[0]]+P[v[1]]+P[v[2]])/4.0;  }};// 增量法求三维凸包。// 注意：没有考虑各种特殊情况（如四点共面）。实践中，请在调用前对输入点进行微小扰动vector<Face> CH3D(Point3* P, int n) {  vector<Face> cur;  // 由于已经进行扰动，前三个点不共线  cur.push_back((Face){{0, 1, 2}});  cur.push_back((Face){{2, 1, 0}});  for(int i = 3; i < n; i++) {    vector<Face> next;    // 计算每条边的“左面”的可见性    for(int j = 0; j < cur.size(); j++) {      Face& f = cur[j];      int res = f.cansee(P, i);      if(!res) next.push_back(f);      for(int k = 0; k < 3; k++) vis[f.v[k]][f.v[(k+1)%3]] = res;    }    for(int j = 0; j < cur.size(); j++)      for(int k = 0; k < 3; k++) {        int a = cur[j].v[k], b = cur[j].v[(k+1)%3];        if(vis[a][b] != vis[b][a] && vis[a][b]) // (a,b)是分界线，左边对P[i]可见          next.push_back((Face){{a, b, i}});      }    cur = next;  }  return cur;}// 点在三角形P0, P1, P2中bool PointInTri(Point3 P, Point3 P0, Point3 P1, Point3 P2) {  double area1 = Area2(P, P0, P1);  double area2 = Area2(P, P1, P2);  double area3 = Area2(P, P2, P0);  return dcmp(area1 + area2 + area3 - Area2(P0, P1, P2)) == 0;}// 三角形P0P1P2是否和线段AB相交bool TriSegIntersection(Point3 P0, Point3 P1, Point3 P2, Point3 A, Point3 B, Point3& P) {  Vector3 n = Cross(P1-P0, P2-P0);  Point3 C = B-A;  if(dcmp(Dot(n, B-A)) == 0) return false; // 线段A-B和平面P0P1P2平行或共面  else { // 平面A和直线P1-P2有惟一交点    double t = Dot(n, P0-A) / Dot(n, B-A);    if(dcmp(t) < 0 || dcmp(t-1) > 0) return false; // 不在线段AB上    P = A + (B-A)*t; // 交点    return PointInTri(P, P0, P1, P2);  }}struct ConvexPolyhedron {  int n;  Point3 P[maxn], P2[maxn];  vector<Face> faces;  ConvexPolyhedron() {    n = 0;  }  void copy_from(const ConvexPolyhedron& other) {    n = other.n;    for(int i = 0; i < n; i++) { P[i] = other.P[i]; P2[i] = other.P2[i]; }    faces = other.faces;  }  void build() {    sort(P, P+n);    n = unique(P, P+n) - P;    for(int i = 0; i < n; i++) P2[i] = add_noise(P[i]);    faces.clear();    if(n >= 4) faces = CH3D(P2, n);  }  double volume() {    if(n < 4) return 0.0;    double res = 0;    for(int i = 0; i < faces.size(); i++)      res += faces[i].Volume6(P, P[0]);    return -res / 6.0;  }  bool is_inside(const Point3& p) {    if(n < 4) return false;    for(int i = 0; i < faces.size(); i++)      if(dcmp(faces[i].Volume6(P, p)) > 0) return false;    return true;  }  void intersect(ConvexPolyhedron& other) {    Point3 inter[maxn];    Point3 T1[3], T2[3], PP;    // get intersections    int c = 0;    for(int f1 = 0; f1 < faces.size(); f1++) {      for(int i = 0; i < 3; i++) T1[i] = P2[faces[f1].v[i]];      for(int f2 = 0; f2 < other.faces.size(); f2++) {        for(int i = 0; i < 3; i++) T2[i] = other.P2[other.faces[f2].v[i]];        for(int i = 0; i < 3; i++) {          if(TriSegIntersection(T1[0], T1[1], T1[2], T2[i], T2[(i+1)%3], PP)) inter[c++] = PP;          if(TriSegIntersection(T2[0], T2[1], T2[2], T1[i], T1[(i+1)%3], PP)) inter[c++] = PP;        }      }    }    for(int i = 0; i < n; i++)      if(other.is_inside(P[i])) inter[c++] = P[i];    for(int i = 0; i < other.n; i++)      if(is_inside(other.P[i])) inter[c++] = other.P[i];    if(c) {      sort(inter, inter+c);      c = unique(inter, inter+c) - inter;    }    // rebuild    n = c;    for(int i = 0; i < n; i++) P[i] = inter[i];    build();  }};#include<cassert>const int maxt = 10;int n, k, m;ConvexPolyhedron CH[maxt];double intersectS(int S) {  ConvexPolyhedron res;  int i;  for(i = 0; i < n; i++) if(S & (1<<i)) {    res.copy_from(CH[i]);    break;  }  for(i++; i < n; i++) if(S & (1<<i))    res.intersect(CH[i]);  return res.volume();}int C(int n, int m) {  int ans = 1;  for(int i = 0; i < m; i++) ans *= n-i;  for(int i = 1; i <= m; i++) ans /= i;  return ans;}double volumeK(int k) {  double ans = 0.0;  for(int S = 0; S < (1<<n); S++) {    int c = 0;    for(int i = 0; i < n; i++) if(S & (1<<i)) c++;    if(c < k) continue;    double vol = intersectS(S);    if((c - k) % 2 == 1) ans -= vol * C(c,k); else ans += vol * C(c,k);  }  return ans;}int main() {  int T;  scanf("%d", &T);  while(T--) {    scanf("%d%d%d", &n, &k, &m);    for(int i = 0; i < n; i++) CH[i].n = 0;    for(int i = 0; i < m; i++) {      int t;      double x, y, z;      scanf("%d%lf%lf%lf", &t, &x, &y, &z); t--;      CH[t].P[CH[t].n++] = Point3(x, y, z);    }    for(int i = 0; i < n; i++) {      int n2 = CH[i].n;      CH[i].build();      assert(n2 == CH[i].n);    }    double ans = 0.0;    for(int i = k; i <= n; i++) ans += volumeK(i);    printf("%.5lf\n", ans);  }  return 0;}