2021牛客暑期多校训练营2-Girlfriend 计算几何

Girlfriend

solution

$|AP|\ge k|BP|,(k>1)$，在二维平面上称为阿波罗尼斯圆，此题三维即阿波罗尼斯球.

设$A(x_0,y_0,z_0),B(x_1,y_1,z_1),P(x,y,z)$，我们有（下面均取等号位置），

$(x-x_0{)}^2+(y-y_0{)}^2+(z-z_0{)}^2=k^2(x-x_1{)}^2+k^2(y-y_1{)}^2+k^2(z-z_1{)}^2$

$=>(x+\frac{x_0-k^2{x}_1}{k^2-1}{)}^2+(y+\frac{y_0-k^2{y}_1}{k^2-1}{)}^2+(z+\frac{z_0-k^2{z}_1}{k^2-1}{)}^2$

$=\frac{(x_0-k^2{x}_1{)}^2+(y_0-k^2{y}_1{)}^2+(z_0-k^2{z}_1{)}^2}{(k^2-1{)}^2}-\frac{(k^2{x}_1^2-x_0^2)+(k^2{y}_1^2-y_0^2)+(k^2{z}_1^2-z_0^2)}{k^2-1}$

同理求得另一个球的方程，再取两球的相交部分即可.

code

#include <bits/stdc++.h>
#define CLR(a,b) memset(a,b,sizeof(a));
const int inf=0x3f3f3f3f;
using namespace std;
const double PI = acos(-1.0);
typedef unsigned long long ll;
const int maxn= 110;
typedef struct point {
    double x,y,z;
    point() {}
    point(double a, double b,double c) {x = a;y = b;z = c;}
    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;
    }
    void output() {
        cout << x << ' ' << y << ' ' << z << '\n';
    }
}point;

double dist(point p1, point p2) { 
    return sqrt((p1 - p2)*(p1 - p2));
}
typedef struct sphere {
    double r;
    point centre;
    void output() {
        centre.output();
        cout << r << '\n';
    }
}sphere;
sphere s,a[10];
void SphereInterVS(sphere a, sphere b,double &v,double &s) {
    double d = dist(a.centre, b.centre);
    double t = (d*d + a.r*a.r - b.r*b.r) / (2.0 * 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.0 * a.r*d));
    double angle_b = 2 * acos((b.r*b.r + d*d - a.r*a.r) / (2.0 * 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;
}
int t, n;
double x, y, z, r;
int cas = 1;
int main() {
//     freopen("input", "r", stdin);
    cin >> t;
    while(t--) {
        double x0, y0, z0, x1, y1, z1;
        double x2, y2, z2, x3, y3, z3;
        double k1, k2;
        cin >> x0 >> y0 >> z0;
        cin >> x1 >> y1 >> z1;
        cin >> x2 >> y2 >> z2;
        cin >> x3 >> y3 >> z3;
        cin >> k1 >> k2;
 
        a[1].centre = {-(x0 - k1 * k1 * x1) / (k1 * k1 - 1),
                       -(y0 - k1 * k1 * y1) / (k1 * k1 - 1),
                       -(z0 - k1 * k1 * z1) / (k1 * k1 - 1)};
        s.centre = {-(x2 - k2 * k2 * x3) / (k2 * k2 - 1),
                       -(y2 - k2 * k2 * y3) / (k2 * k2 - 1),
                       -(z2 - k2 * k2 * z3) / (k2 * k2 - 1)};
 
 
 
        a[1].r = ((k1 * k1 * x1 * x1 - x0 * x0) +
               (k1 * k1 * y1 * y1 - y0 * y0) +
               (k1 * k1 * z1 * z1 - z0 * z0)) / (k1 * k1 - 1) -
              ((x0 - k1 * k1 * x1) * (x0 - k1 * k1 * x1) +
               (y0 - k1 * k1 * y1) * (y0 - k1 * k1 * y1) +
               (z0 - k1 * k1 * z1) * (z0 - k1 * k1 * z1)) / ((k1 * k1 - 1) * (k1 * k1 - 1));
 
        s.r = ((k2 * k2 * x3 * x3 - x2 * x2) +
               (k2 * k2 * y3 * y3 - y2 * y2) +
               (k2 * k2 * z3 * z3 - z2 * z2)) / (k2 * k2 - 1) -
              ((x2 - k2 * k2 * x3) * (x2 - k2 * k2 * x3) +
               (y2 - k2 * k2 * y3) * (y2 - k2 * k2 * y3) +
               (z2 - k2 * k2 * z3) * (z2 - k2 * k2 * z3)) / ((k2 * k2 - 1) * (k2 * k2 - 1));
 
        a[1].r = sqrt(-a[1].r);
        s.r = sqrt(-s.r);
 
//        a[1].output();
//        s.output();
 
        double ans = 0, v = 0;
 
        for(int i = 1; i <= 1; i++) {
            double ss, dis = dist(s.centre, a[i].centre);
 
            if(dis >= s.r + a[i].r)continue;  //在外部
            if(dis + min(s.r, a[i].r) <= max(s.r, a[i].r))  //在内部
            {
                ans += 4.0 / 3.0 * PI * min(s.r,a[i].r) * min(s.r,a[i].r) * min(s.r,a[i].r);
                continue;
            }
            SphereInterVS(s, a[i], v, ss); //相交部分
            ans += v;
        }
        cout << fixed << setprecision(3) << ans << '\n';
    }
}