湖南省第八届程序设计竞赛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;
}