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}^noperatorname{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}。
1le nle 20000,1le mle5,0le 所有坐标 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_ivec Q_i。
由于 text{dist}(A,B)=sqrt{(A-B)^2}(这里的平方指自己与自己做点乘),因此可以对于每个 kin[1,t) 列出方程:
拆平方移个项即可得到:
用高斯消元解个方程即可。
Code
代码语言:javascript复制#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;
}