2021牛客暑期多校训练营2 Girlfriend（阿波罗尼斯圆 + 计算几何）

题意

空间里有6个点，满足$\left|P_{1}A\right|⩾k_1\left|P_{1}B\right|$，$\left|P_{2}C\right|⩾k_2\left|P_{2}D\right|$，求P₁，P₂各自的轨迹相交部分的体积；

思路

对于固定的k，P1，P2的轨迹构成标准的阿波罗尼斯球壳；
同时，题意的k > 1，可以得知为轨迹为实心的球体，于是问题转化为求两球相交部分的体积；

将条件$\left|P_{1}A\right|⩾k_1\left|P_{1}B\right|$，$\left|P_{2}C\right|⩾k_2\left|P_{2}D\right|$化简得：

$$\left(1-k^2\right)\left(x_p+y_p^2+z_p^2\right)+x_p\left(2k^2{x}_2-2x_1\right)+y_p\left(2k^2{y}_2-2y_1\right)+z_p\left(2k^2{z}_2-2z_1\right)+x_1^2+y_1^2+z_1^2-k^2\left(x_2^2+y_2^2+z_2^2\right)=0$$

又知球的方程：$x^2+y^2+z^2+2ax+2by+2cz+d=0$

可得球心坐标$\left(-a,-b,-c\right)$，半径$r=\sqrt{a^2+b^2+c^2-d}$

然后判断两个球是否相交即可；

如果相交

把立体具化成平面，黄色是相交部分在球1中的体积，利用三重积分或者二重积分算体积；
$$\int \int \int dxdydz={\int }_{r-h}^r\pi \ast (r^2-z^2)dz$$

或者

$$\int \int \left[\sqrt{r^2-x^2-y^2}-(r-h)\right]dxdy={\int }_0^{2\pi }d\Theta {\int }_0^{\sqrt{2rh-h^2}}\left(\sqrt{r^2-{\rho }^2}-\left(r-h\right)\right)\rho d\rho$$

化简可得到
$$V=\pi h^2\left(r-\frac{1}{3}h\right)$$

最后同理得到球2部分加起来即可；

AC代码

#include <iostream>
#include <cstring>
#include <cmath>
#include <cstdio>
#include <string>
#include <algorithm>
#include <queue>
#include <utility>
#include <stack>
#include <map>
#include <vector>
#include <set>
#include <iomanip>
#include <unordered_map>
#define hz020 return
#define mes memset
#define mec memcpy
#define x first
#define y second

using namespace std;

typedef long long ll;
typedef unsigned long long ull;
typedef pair<int, int> pii;
typedef pair<ll,ll> PII;

const double pi = acos(-1);
const int N = 25;
const int null = 0x3f3f3f3f,INF = 1e9;
const ll mod = 998244353;

int T;
double x[4], y[4], z[4];
int k1, k2;

void solve(double x1, double y1, double z1, double x2, double y2, double z2, double r1, double r2)
{
	double dist = sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) + (z1 - z2) * (z1 - z2));
	
	if (dist >= r1 + r2) printf("0.000\n");
	else if (dist + r2 <= r1) printf("%.3lf\n", pi * r2 * r2 * r2 * 4 / 3);
	else if (dist + r1 <= r2) printf("%.3lf\n", pi * r1 * r1 * r1 * 4 / 3);
	else
	{
		double cos1 = (dist * dist + r1 * r1 - r2 * r2) / (2 * r1 * dist), cos2 = (dist * dist + r2 * r2 - r1 * r1) / (2 * r2 * dist);
		double h1 = r1 - r1 * cos1, h2 = r2 - r2 * cos2;
		
		printf("%.3lf\n", pi * h1 * h1 * (r1 - h1 / 3) + pi * h2 * h2 * (r2 - h2 / 3));
	}
}

int main()
{
	cin >> T;
	
	while (T --)
	{
		for (int i = 0; i <= 3;i ++ ) cin >> x[i] >> y[i] >> z[i];
		cin >> k1 >> k2;
		
		double sq1 = k1 * k1, idx1 = 1 - sq1;
		//圆心 
		double xp1 = (sq1 * x[1] - x[0]) / idx1, yp1 = (sq1 * y[1] - y[0]) / idx1, zp1 = (sq1 * z[1] - z[0]) / idx1;
	
		double dp1 = (x[0] * x[0] + y[0] * y[0] + z[0] * z[0] - sq1 * (x[1] * x[1] + y[1] * y[1] + z[1] * z[1])) / idx1;
		//半径 
		double rp1 = sqrt(xp1 * xp1 + yp1 * yp1 + zp1 * zp1 - dp1);
		
		double sq2 = k2 * k2, idx2 = 1 - sq2;
		//圆心 
		double xp2 = (sq2 * x[3] - x[2]) / idx2, yp2 = (sq2 * y[3] - y[2]) / idx2, zp2 = (sq2 * z[3] - z[2]) / idx2;
	
		double dp2 = (x[2] * x[2] + y[2] * y[2] + z[2] * z[2] - sq2 * (x[3] * x[3] + y[3] * y[3] + z[3] * z[3])) / idx2;
		//半径 
		double rp2 = sqrt(xp2 * xp2 + yp2 * yp2 + zp2 * zp2 - dp2);
		
		solve (xp1, yp1, zp1, xp2, yp2, zp2, rp1, rp2);
	}

	
	hz020 0;
}  

END