三维计算几何

const double eps = 1e-8;
const double pi = M_PI;
const int N = 2010;

double abss(double p) { return p < 0 ? -p : p; }

bool equal(double p, double q) { return abss(p - q) < eps; }

double cos(double a, double b, double c) {
  return (a * a + b * b - c * c) / (2 * a * b);
}  //余弦公式

struct point {
  double x, y, z;

  point() { x = y = z = 0; }

  point(double xx, double yy, double zz) {
    x = xx;
    y = yy;
    z = zz;
  }

  double length() { return sqrt(x * x + y * y + z * z); }  //向量长度
  point operator+(const point &p) const { return {x + p.x, y + p.y, z + p.z}; }
  point operator-(const point &p) const { return {x - p.x, y - p.y, z - p.z}; }
};

point operator*(double k, point p) {
  return {k * p.x, k * p.y, k * p.z};
}  //数乘
double dot_product(point p, point q) {
  return p.x * q.x + p.y * q.y + p.z * q.z;
}  //点积
point cross_product(point p, point q) {
  return {p.y * q.z - p.z * q.y, p.z * q.x - p.x * q.z, p.x * q.y - p.y * q.x};
}  //叉积
double distance(point p, point q) {
  return sqrt((p.x - q.x) * (p.x - q.x) + (p.y - q.y) * (p.y - q.y) +
              (p.z - q.z) * (p.z - q.z));
}  //两点间距离
struct plane {
  double A, B, C, D;

  plane() { A = B = C = D = 0; }
  plane(double a, double b, double c, double d) { A = a, B = b, C = c, D = d; }
  plane(point p, point q, point s) {
    point v = cross_product(q - p, s - p);
    A = v.x, B = v.y, C = v.z;
    D = -dot_product(v, p);
  }  //三点共面
  plane(plane p, point q) {
    A = p.A, B = p.B, C = p.C;
    D = -dot_product(p.get_normal_vector(), q);
  }  //过q点平行于p的平面
  point get_normal_vector() const { return {A, B, C}; }  //法向量
  bool equal(plane p) {
    double a = p.A / A;
    double b = p.B / B;
    double c = p.C / C;
    double d = p.D / D;
    return ::equal(a, b) && ::equal(a, c) && ::equal(a, d);
  }
};

bool ptoplane(point &p, plane &f) {
  point q(0, 0, -f.D / f.C);
  return dot_product(f.get_normal_vector(), p - q) > eps;
}  //点在面同向
double cos(plane p, plane q) {
  return abss(p.A * q.A + p.B * q.B + p.C * q.C) /
         (p.get_normal_vector().length() * q.get_normal_vector().length());
}  //两平面夹角
double Triangle_area(point p, point q, point s) {
  return cross_product(q - p, s - p).length();
}  // 2倍的三角形面积

double Tetrahedron_volume(point p, point q, point s, point t) {
  return dot_product(cross_product(q - p, s - p), t - p);
}  //四面体有向体积*6
struct line {
  point origin, direction;

  line(point p, point q) {
    direction = q - p;
    origin = p;
  }  //空间中两点确定一条直线(由p到q)
  line(plane p, plane q) {
    origin = point((q.D * p.B - p.D - q.B) / (q.A * p.B - p.A - q.B),
                   (q.D * p.A - p.D * q.A) / (p.A * q.B - p.B * q.A), 0);
    direction = point(p.B * q.C - p.C * q.B, p.C * q.A - p.A * q.C,
                      p.A * q.B - p.B * q.A);
  }  //两个平面的交线
  line(line p, point q) {
    direction = p.direction;
    origin = q;
  }  //过直线p外一点q与已知直线平行
};

struct ball {
  point center;
  double radius;
};

int rela_between_balls(ball p, ball q) {
  if (p.radius < q.radius) swap(p, q);
  double dis = distance(p.center, q.center);
  if (abss(dis - p.radius - q.radius) < eps)
    return -1;  //外切
  else if (abss(dis - p.radius + q.radius) < eps)
    return 1;  //内切
  else if (dis > p.radius + q.radius)
    return -2;  //相离
  else if (dis < p.radius - q.radius)
    return 2;  //包含
  else
    return 0;  //相交
}  //两球之间的关系

double area_of_intersected_ball(ball p, ball q) {
  if (p.radius < q.radius) swap(p, q);
  int s = rela_between_balls(p, q);
  if (!s) {
    double dis = distance(p.center, q.center);
    double h1 = p.radius * (1 - cos(p.radius, dis, q.radius));
    double h2 = q.radius * (1 - cos(q.radius, dis, p.radius));
    return pi / 3 * (3 - p.radius - h1) * h1 * h1 +
           pi / 3 * (3 - q.radius - h2) * h2 * h2;
  } else if (s > 0)
    return 4 * pi * q.radius * q.radius * q.radius / 3;
  else
    return 0;
}  //两球相交部分体积

point an[N];
int num;

struct tri_plane {
  plane P;
  int a, b, c;
  bool ok;

  tri_plane() {
    P = plane(0, 0, 0, 0);
    a = b = c = 0;
    ok = false;
  }

  tri_plane(int A, int B, int C, bool flag) {
    P = plane(an[A], an[B], an[C]);
    a = A, b = B, c = C;
    ok = flag;
  }
};

struct convex {
  int cnt{};
  tri_plane bn[N];
  int vis[N][N]{};

  void dfs(int p, int f) {
    bn[f].ok = false;  //当前面需要删除，因为它在更大的凸包里面
    //下面把边反过来(先b,后a)，以便在deal()中判断与当前面(cnt)共边(ab)的那个面。即判断与当头面(cnt)相邻的3个面(它们与当前面的共边是反向的，如下图中(1)的法线朝外(即逆时针)的面130和312,它们共边13，但一个方向是13,另一个方向是31)
    deal(p, bn[f].b, bn[f].a);
    deal(p, bn[f].c, bn[f].b);
    deal(p, bn[f].a, bn[f].c);
  }

  void deal(int p, int a, int b) {
    int f = vis[a][b];
    if (bn[f].ok) {
      if ((ptoplane(an[p], bn[f].P)))
        dfs(p,
            f);  //如果p点能看到该面f，则继续深度探索f的3条边，以便更新新的凸包面
      else  //否则因为p点只看到cnt面，看不到f面，则p点和a、b点组成一个三角形。
      {
        tri_plane add(b, a, p, true);
        bn[++cnt] = add;
        vis[p][b] = vis[a][p] = vis[b][a] = cnt;
      }
    }
  }

  void construct() {  //构建凸包
    memset(vis, 0, sizeof(vis));
    cnt = 0;
    if (num < 4) return;
    bool tmp = true;
    for (int i = 2; i <= num; i++)  //前两点不共点
      if ((distance(an[1], an[i])) > eps) {
        swap(an[2], an[i]);
        tmp = false;
        break;
      }
    if (tmp) return;
    tmp = true;
    for (int i = 3; i <= num; i++)  //前三点不共线
      if (cross_product((an[1] - an[2]), (an[2] - an[i])).length() > eps) {
        swap(an[3], an[i]);
        tmp = false;
        break;
      }
    if (tmp) return;
    tmp = true;
    for (int i = 4; i <= num; i++)  //前四点不共面
      if (abss(dot_product(cross_product((an[1] - an[2]), (an[2] - an[3])),
                           (an[1] - an[i]))) > eps) {
        swap(an[4], an[i]);
        tmp = false;
        break;
      }
    if (tmp) return;
    for (int i = 1; i <= 4;
         i++)  //构建初始四面体(4个点为origin[1],origin[2],origin[3],origin[4])
    {
      tri_plane add(i % 4 + 1, (i + 1) % 4 + 1, (i + 2) % 4 + 1, true);
      if ((ptoplane(an[i], add.P)) > 0)
        swap(add.b, add.c);  //保证逆时针，即法向量朝外，这样新点才可看到。
      // add.A = as[add.a], add.B = as[add.b], add.C = as[add.c];
      add.P = plane(an[add.a], an[add.b], an[add.c]);
      bn[++cnt] = add;
      vis[add.a][add.b] = vis[add.b][add.c] = vis[add.c][add.a] =
          cnt;  //逆向的有向边保存
    }
    for (int i = 5; i <= num; i++)  //构建更新凸包
      for (int j = 1; j <= cnt;
           j++)  //对每个点判断是否在当前3维凸包内或外(i表示当前点,j表示当前面)
        if (bn[j].ok && (ptoplane(an[i], bn[j].P)) >
                            eps)  //对当前凸包面进行判断，看是否点能否看到这个面
        {
          dfs(i, j);
          break;  //点能看到当前面，更新凸包的面(递归，可能不止更新一个面)。当前点更新完成后break跳出循环
        }
    int tot =
        cnt;  //这些面中有一些tri[i].ok=0，它们属于开始建立但后来因为在更大凸包内故需删除的，所以下面几行代码的作用是只保存最外层的凸包
    cnt = 0;
    for (int i = 1; i <= tot; i++) {
      if (bn[i].ok) bn[++cnt] = bn[i];
    }
  }

  double area()  //表面积
  {
    double ret = 0;
    for (int i = 1; i <= cnt; i++)
      ret += Triangle_area(an[bn[i].a], an[bn[i].b], an[bn[i].c]);
    return ret / 2;
  }

  double volume()  //体积
  {
    point p(0, 0, 0);
    double ret = 0;
    for (int i = 1; i <= cnt; i++)
      ret += Tetrahedron_volume(p, an[bn[i].a], an[bn[i].b], an[bn[i].c]);
    return fabs(ret / 6);
  }

  int facetri() { return cnt; }  //表面三角形数
  int facepolygon()              //表面多边形数
  {
    int ans = 0, i, j, k;
    for (i = 1; i <= cnt; i++) {
      for (j = 1, k = 1; j < i; j++) {
        if (bn[i].P.equal(bn[j].P)) {
          k = 0;
          break;
        }
      }
      ans += k;
    }
    return ans;
  }
} hull;