牛客2021暑假多校训练(2) F - 计算几何 两球之交 + 定比分点 + 阿波罗尼斯圆性质

题意

• 组询问，每次给定空间中的四个定点，以及两个常数，动点制约于.求构成的几何体的体积交，答案精度误差不得超过.

  制约

• ，保证互异。

  思路

  高中数学复习：对于两定点，动点满足，则有： 是一条直线 是 一个圆，称作阿波罗尼斯圆，简称阿氏圆

  证明：设.以中点为坐标原点，作为轴建立平面直角坐标系，则.设P,由,有，两边平方，化简得
  是的中垂线

• 而在三维空间中，其形态自然是一个球。这实际上是阿氏圆推出来的，后面的两个球方程该怎么搞呢？按照上面的建系的方式来推，避免不了要化简复杂的算式，所以下面要介绍一个公式(实际上也是数学手的第一直觉)：定比分点公式

  

  多么优美的算式

由此，我们就可以确定所在点的最小值了(取等号时刻的位置)，那么最远处该怎么算呢？我门取与之相对称的位置(负方向的等号)即可：。圆心在两点的中点处，半径则是两点距离的一半，剩下只需要判断相离(相切)、内含(外含)与相交的情况了套板子。

模板+代码

• 模板推导过程

#include <bits/stdc++.h>
using namespace std;
static const double PI = acos(-1.0);
struct point {
    double x,y,z;
    point(double _x = .0, double _y = .0,double _z = .0) 
    : x(_x), y(_y), z(_z) { }
    point operator -(const point &b)const { return point(x - b.x, y - b.y, z - b.z); }
    point operator +(const point &b)const { return point(x + b.x, y + b.y, z + b.z); }
    point operator *(const double &k)const { return point(x * k, y * k,z * k); }
    point operator /(const double &k)const { return point(x / k, y / k, z / k); }
    double operator *(const point &b)const { return x * b.x + y * b.y + z * b.z; }
    friend point operator*(double x, point p) { return p * x; }
};

double dist(point p1, point p2) {return sqrt((p1 - p2) * (p1 - p2)); }

struct sphere {
    point center;
    double r;
    sphere(point _center = point(), double _r = 1.)
    : center(_center), r(_r) { }
};

void SphereInterVS(sphere a, sphere b,double &v,double &s) {
    double d = dist(a.center, b.center);
    double t = (d * d + a.r * a.r - b.r * b.r) / (2. * d);
    double h = sqrt((a.r * a.r) - (t * t)) * 2;
    double angle_a = 2 * acos((a.r*a.r + d*d - b.r*b.r) / (2. * a.r*d));
    double angle_b = 2 * acos((b.r*b.r + d*d - a.r*a.r) / (2. * b.r*d));
    double l1 = ((a.r * a.r - b.r * b.r) / d + d) / 2;
    double l2 = d - l1;
    double x1 = a.r - l1, x2 = b.r - l2;
    double v1 = PI * x1 * x1 * (a.r - x1 / 3);
    double v2 = PI * x2 * x2 * (b.r - x2 / 3);
    v = v1 + v2;
    double s1 = PI * a.r * x1;
    double s2 = PI * a.r * x2;
    s = 4. * PI * (a.r*a.r + b.r*b.r) - s1 - s2;
}

sphere set_sphere(point a, point b, double k) {
    point t1(1. / (1 + k) * a + 1. * k / (1 + k) * b);
    point t2(1. * k / (k - 1) * b - 1. / (k - 1) * a);
    return sphere((t1 + t2) / 2., dist(t1, t2) / 2.);
}
int main() {
    // ios::sync_with_stdio(false), cin.tie(nullptr);
    int tt; cin >> tt; while (tt --) {
        sphere a, b;
        int k1, k2, x1, x2, y1, y2, z1, z2, x3, x4, y3, y4, z3, z4;
        double v = 0., s = .0;

        cin >> x1 >> y1 >> z1; point p1(x1, y1, z1);
        cin >> x2 >> y2 >> z2; point p2(x2, y2, z2);
        cin >> x3 >> y3 >> z3; point p3(x3, y3, z3);
        cin >> x4 >> y4 >> z4; point p4(x4, y4, z4);

        cin >> k1 >> k2;

        a = set_sphere(p1, p2, k1);
        b = set_sphere(p3, p4, k2);

        double ans = .0, dis = dist(a.center, b.center);

        if(dis >= a.r + b.r) { // 两球相离
            ans = .0;
        } else if(dis + min(a.r, b.r) <= max(a.r, b.r)) { // 两球内含
            ans = 4. / 3. * PI * min(a.r, b.r) * min(a.r, b.r) * min(a.r, b.r);
        } else {
            SphereInterVS(a, b, v, s);
            ans = v;
        }
        printf("%.10f\n", ans);
    }
    return 0;
}