YbtOJ 894「高斯消元」高维寻点
YbtOJ 894「高斯消元」高维寻点
题目链接：YbtOJ #894
小 A 有一个 m 维空间。在这个空间中有 n 个特殊点，其中第 i 个特殊点 p_i 的坐标为 (x_{i,1},x_{i,2},\cdots,x_{i,m})。
他希望找到这个 m 维空间中的一个点 P，使得 \max_{i=1}^n\operatorname{dist}(P,p_{i}) 最小。
提示：m 维空间中的两点 (A_1,A_2,\cdots,A_m),(B_1,B_2,\cdots,B_m) 间的距离为 \sqrt{\sum_{i=1}^m(A_i-B_i)^2}。
1\le n\le 20000，1\le m\le5，0\le 所有坐标 \le10^4。
Solution
首先二维最小覆盖圆的求法：
首先我们枚举一个点 p_i，如果它不在原本前 i-1 个点的最小覆盖圆内，就必然在当前前 i 个点的最小覆盖圆上。因此我们重构最小覆盖圆，由于初始只能确定这一个点在最小覆盖圆上，所以令此时的最小覆盖圆的圆心为当前点，半径为0。
然后我们再在 [1,i) 中枚举点 p_j，如果它不在当前的最小覆盖圆内，同理我们可以确定它在新的最小覆盖圆上，那么我们就能确定两个点。所以令此时的最小覆盖圆的圆心为这两个点构成线段的中点，半径就是这两点间距离的一半。
同理继续在 [1,j) 中枚举点 p_k，如果它不在当前的最小覆盖圆内，就令新的最小覆盖圆为这三个点的最小外接圆（其实之前两种情况也都是特殊的最小外接圆）。
这个做法看似暴力，但实际上可以利用 随机增量法 来使复杂度正确。即将数据随机打乱。
可以证明是 O(N) 的。
那么高维的只需要解决如何求最小外接圆。
令 \vec Q_i=q_i-q_t，设圆心 O=q_t+\sum_{i=1}^{t-1}\lambda_i\vec Q_i。
由于 \text{dist}(A,B)=\sqrt{(A-B)^2}（这里的平方指自己与自己做点乘），因此可以对于每个 k\in[1,t) 列出方程：
拆平方移个项即可得到：
用高斯消元解个方程即可。
Code
#pragma GCC optimize("Ofast")
#pragma GCC target("sse,sse2,sse3,ssse3,sse4,popcnt,abm,mmx,avx,avx2,fma")
#pragma GCC optimize("unroll-loops")
#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define W while
#define I inline
#define RI register int
#define LL long long
#define Cn const
#define CI Cn int&
using namespace std;
namespace Debug{
Tp I void _debug(Cn char* f,Ty t){cerr<<f<<'='<<t<<endl;}
Ts I void _debug(Cn char* f,Ty x,Ar... y){W(*f!=',') cerr<<*f++;cerr<<'='<<x<<",";_debug(f+1,y...);}
Tp ostream& operator<<(ostream& os,Cn vector<Ty>& V){os<<"[";for(Cn auto& vv:V) os<<vv<<",";os<<"]";return os;}
#define gdb(...) _debug(#__VA_ARGS__,__VA_ARGS__)
}using namespace Debug;
Cn int N=2e4+10,M=6;
int n,m;
double A[M][M],b[M];
#define P valarray<double>
P a[N];
vector<P> v,t;
#define pb push_back
struct C{P O;double R;}Ans;
I bool IC(Cn C& x,Cn P& p){return ((x.O-p)*(x.O-p)).sum()<=x.R;}
C G(){
RI i,j,k,sz=v.size();for(t=v,memset(A,0,sizeof(A)),i=0;i<sz-1;i++) t[i]-=t.back();for(i=0;i<sz-1;i++) for(A[i][sz-1]=(t[i]*t[i]).sum(),j=0;j<sz-1;j++) A[i][j]=2*(t[i]*t[j]).sum();
for(i=0;i<sz-1;i++){for(j=i;j<sz-1;j++) if(fabs(A[j][i])>fabs(A[i][i])) swap(A[i],A[j]);for(j=sz-1;j>=i;j--) A[i][j]/=A[i][i];for(j=0;j<sz-1;j++) if(i^j){double t=A[j][i]/A[i][i];for(k=0;k<sz;k++) A[j][k]-=A[i][k]*t;}}
P S;S.resize(m);for(i=0;i<sz-1;i++) S=S+A[i][sz-1]*t[i];return (C){t[sz-1]+S,(S*S).sum()};
}
I C Sol(CI x){C t;t.O.resize(m),t.R=0,v.size()&&(t=G(),0);RI i;for(i=1;i<=x;i++) !IC(t,a[i])&&(v.pb(a[i]),t=Sol(i-1),v.pop_back(),0);return t;}
int main(){
freopen("dimension.in","r",stdin),freopen("dimension.out","w",stdout);
RI i,j;double x;for(scanf("%d%d",&n,&m),i=1;i<=n;i++) for(a[i]=P(m),j=0;j<m;j++) scanf("%lf",&x),a[i][j]=x;
for(Ans=Sol(n),i=0;i<m;i++) printf("%.10lf ",Ans.O[i]);return printf("\n"),0;
}