【计算几何，数学】7.14 T3 @ xdfz

Problem Link

给定 $n$ 个球和一个点 $P$，求点 $P$ 到这些球的交内一点的距离的最小值。保证有解。$n\le 10^6$。

---

和最小圆覆盖一个套路。考虑维护一个当前答案，初始即为询问点 $P$。

从 $1$ 到 $n$ 枚举 $i$，如果当前答案在球 $i$ 内则跳过，否则可以证明答案一定在 $i$ 的球面上。证明后面再说。

做法是：直 remake，将当前答案设为 $P$ 到球 $i$ 表面的最小距离点，然后枚举 $j=1\sim i-1$，和上面一样，如果当前答案在球 $j$ 内，那无所谓；否则，同理可以说明答案也一定在 $j$ 的球面上（注意此时答案在 $i$ 的球面上的前提仍然保持），所以推出答案在球 $i$ 和球 $j$ 相交所得的圆上，找出 $P$ 到这个圆的距离作为最新答案。然后继续枚举 $k=1\sim j-1$，同理如果当前答案不在球 $k$ 内则更新为 $i,j,k$ 三个球的两个交点之中离 $P$ 较近的那个。

首先来考虑正确性。当我们枚举 $i$ 发现出问题的时候，由于球交是凸的，所以原先答案到此时的答案（此时的答案定义为 $1\sim i$ 中某 3 个球的交点的最优解）一定存在一条连续的移动路径使得移动过程中和 $P$ 的距离始终不变，于是最优解一定在 $i$ 的球面上（而不是球内）。

然后枚举 $j$，当前答案定义为“考虑 $i$ 和 $\le j$ 的某两个球交出的点的最优答案”，那么如果 $j$ 不包含最优解，同上是可以移动到 $j$ 的球面上的。$k$ 同理。

至于为什么答案一定可以由 $1\sim 3$ 个球的球交中的最优解得到，可以考虑最终答案一定在某个“犄角旮旯”，那里是由至多 3 个球交出来的。

那么时间复杂度呢？和最小圆覆盖一样分析，开始的时候随机重排所有点，

对于单个 $j$，枚举 $k$ 的时间复杂度是 $O(j)$ 的。对于单个 $i$，考虑每个 $1\le j<i$，$j$ 能更新答案，意味着 $j$ 是 “$i$ 和 $1\sim j$ 中某两个球能组成答案的最优解”所取到的那两个球之一，这个概率是 $2/j$！所对单个 $i$ 时间复杂度是 $\sum_{j=1}^{i-1}O(j)\cdot 2/j =O(i)$。用完全相同的分析可以得到总时间复杂度为 $O(n)$。

点击查看代码

#include <bits/stdc++.h>
#define For(i,a,b) for(int i=a;i<=b;i++)
#define Rev(i,a,b) for(int i=a;i>=b;i--)
#define Fin(file) freopen(file,"r",stdin);
#define Fout(file) freopen(file,"w",stdout);
using namespace std;
const int N=1e6+5; using ll = long long; const double EPS=1e-6;
struct Point{
	double x,y,z;
	Point(double xx=0,double yy=0,double zz=0) : x(xx),y(yy),z(zz) {}
};
typedef Point Vector;
inline Vector operator- (Point A,Point B) { return Vector(A.x-B.x,A.y-B.y,A.z-B.z); }
inline Point operator+ (Point A,Vector v) { return Point(A.x+v.x,A.y+v.y,A.z+v.z); }
inline Vector operator* (Vector v,double o) { return Vector(v.x*o,v.y*o,v.z*o); }
inline Vector operator/ (Vector v,double o) { return Vector(v.x/o,v.y/o,v.z/o); }
inline double Length2(Vector v) { return v.x*v.x+v.y*v.y+v.z*v.z; }
inline double Dist(Point U,Point V) { return sqrt(Length2(U-V)); }
inline double Dot(Vector a,Vector b) { return a.x*b.x+a.y*b.y+a.z*b.z; }
inline Vector Cross(Vector a,Vector b) { return Vector(a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x); }
inline Point Go(Point A,Point B,double d) { return A+(B-A)*d/Dist(A,B); }
inline int dcmp(double d) { return d<-EPS?-1:d>EPS?1:0; }
int n; Point O[N]; double R[N]; mt19937 rng(190345);
inline bool In(int u,Point A) { return Dist(O[u],A)<=R[u]; }
inline double Cos(double a,double b,double c) { return (a*a+b*b-c*c)/(2*a*b); }
inline Point Perp(Point T,Vector v,Point X){
	auto [A,B,C]=X-T; auto [a,b,c]=v;
	double p=-(a*A+b*B+c*C)/(a*a+b*b+c*c);
	return X+v*p;
}
inline Point solve1(Point U,double R1,Point P) { return Go(U,P,R1); }
inline Point solve2(Point U,Point V,double R1,double R2,Point P){
	double d=Dist(U,V);
	Point T=Go(U,V,R1*Cos(d,R1,R2));
	Point S=Perp(T,V-U,P);
	Point K=Go(T,S,R1*sin(acos(Cos(d,R1,R2))));
	return K;
}
inline Point solve3(Point U,Point V,Point W,double R1,double R2,double R3,Point P){
	double d=Dist(U,V);
	Point T=Go(U,V,R1*Cos(d,R1,R2));
	Point S=Perp(T,V-U,W);
	Point K=Go(T,S,R1*sin(acos(Cos(d,R1,R2))));
	double r1=Dist(T,K),r2=sqrt(R3*R3-Length2(W-S));
	Point X=Go(T,S,r1*Cos(Dist(S,T),r1,r2));
	Vector v=Cross(V-U,W-U); v=v/sqrt(Length2(v));
	double h=sqrt(r1*r1-Length2(T-X));
	Point Y=X+v*h,Z=X-v*h;
	return Dist(Y,P)<Dist(Z,P)?Y:Z;
}
int main(){
	//~ Fin("ball.in"); Fout("ball.out");
	ios::sync_with_stdio(0); cin.tie(0);
	cout<<setprecision(6)<<fixed;
	Point P; cin>>n>>P.x>>P.y>>P.z; For(i,1,n) cin>>O[i].x>>O[i].y>>O[i].z>>R[i];
	For(i,2,n) { int j=rng()%i+1; swap(O[i],O[j]); swap(R[i],R[j]); }
	Point Ans=P;
	For(i,1,n) if(!In(i,Ans)) {
		Ans=solve1(O[i],R[i],P);
		For(j,1,i-1) if(!In(j,Ans)){
			Ans=solve2(O[i],O[j],R[i],R[j],P);
			For(k,1,j-1) if(!In(k,Ans)) {
				Ans=solve3(O[i],O[j],O[k],R[i],R[j],R[k],P);
			}
		}
	}
	cout<<Ans.x<<' '<<Ans.y<<' '<<Ans.z<<'\n';
	return 0;
}