YbtOJ 894「高斯消元」高维寻点

2022-09-19 13:41:49 浏览数 (1)

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 200001le mle50le 所有坐标 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) 列出方程:

(sum_{i=1}^{t-1}lambda_ivec Q_i)^2=(sum_{i=1}^{t-1}lambda_ivec Q_i-vec Q_k)^2

拆平方移个项即可得到:

sum_{i=1}^{t-1}2(vec Q_icdot vec Q_k)lambda_i=(vec Q_k)^2

用高斯消元解个方程即可。

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;
}

0 人点赞