[ICPC2020 WF] What's Our Vector, Victor?

给出 \(d\) 维空间的 \(n\) 个点及它们到某个定点的距离，你需要解出这个定点的坐标。

保证有解，任意输出一组解即可。

\(n,d\le 500\)，所有点的坐标（包括定点）随机。

拜谢 zhy 大师！给出了一个十分简单好理解的做法！线代不好怎么办？我们可以猜猜结论！

首先这道题相当于是给了 \(n\) 个形如 \(\sum_{j=1}^d (a_{i,j}-x_j)^2=b_i^2\) 的方程，我们设 \(\sum_{j=1}^d x_j^2=x_{d+1}\)，那么这道题的方程可以整理成关于 \(d+1\) 维向量 \(\vec{x}\) 的线性方程。

注意到这道题数据纯随机，这个条件在线代题中往往意味着每一个方程都是有用的！即，如果 \(n>d\)，那么我们随便取出其中 \(d+1\) 个方程，组成的矩阵极有可能满秩，解之即可得到 \(\vec{x}\)。

如果 \(n=d\)，难处在于 \(x_{d+1}\) 解不出来，但是我们还有二次方程 \(\sum_{j=1}^d x_j^2=x_{d+1}\)。所以我们可以先把所有 \(x_i(i\le d)\) 表示成关于 \(x_{d+1}\) 的线性表达式，然后带入二次方程中解个一元二次方程。

这道题真正的难点便落在 \(n<d\) 的情况上了。这种情况下，我们发现 \(n\) 个随机 \(d\) 维球壳相交（题目保证有交）结果极大概率是一个 \(d-n+1\) 维的球壳，也就是说我们有无穷多个解，我们需要任取一个，怎么办呢？

考虑对这个方程组进行基变换减少变量的个数。我们找一个一定包含至少一个答案的子线性空间，即选取 \(k\) 个合适的 \(d\) 维向量 \(\vec{v_i}\) 和一个原点 \(\vec{s}\)，将 \(\vec{x}\) 表示成 \(\vec{s}+\sum_{i=1}^k r_i\vec{v_i}\)，那么原线性方程就可以整理成关于 \(k\) 维向量 \(\vec{r}\) 的线性方程。

我们断言，我们可以简单构造一个 \(k=n\) 维的线性子空间一定包含至少两个答案！构造的方法意外的简单，我们先找到包含题目给出的 \(n\) 个球心的极小子线性空间 \(V'\)，由于数据随机 \(V'\) 的维度也极大概率就是 \(n-1\) 维的！然后我们任意寻找一个不在 \(V'\) 里的一个点，考虑同时包含了这个点和 \(V'\) 的极小线性子空间 \(V\)，这个 \(V\) 一定包含了至少两个答案。

严谨证明似乎我和 zhy 都不太会。但是很好进行感性理解，比如说 \(n=1,d=3\) 的情况，\(V'\) 仅包含球心，在 \(V'\) 外任选一个点，\(V\) 就会是一条过圆心的直线，这条直线一定会穿过球壳，交于两个点上。

这个题中如何具体寻找 \(V\) 的基向量 \(\vec{v_i}\) 非常简单！因为随机数据下原坐标系的原点大概率不在 \(V'\) 中，所以取 \(\vec{s}=\vec{0}\)，取 \(\vec{v_i}\) 为原点指向这 \(n\) 个点的向量即可。

有了一个 \(k=n\) 维且包含答案的子空间后，我们直接运行类似 \(n=d\) 情况的算法解二次方程解出 \(\vec{r}\)，就可以还原出 \(\vec{x}\)。

解二次方程时，大概率有两个不重的根，选取哪个根都是可以的！因为 \(V\) 与答案集合大概率恰好有两个交点！

代码十分好写。由于数据随机，高斯消元时你甚至不用移动方程的位置，直接认为第 \(i\) 个方程主元是 \(i\) 就可以了！

注意这份代码无法通过样例，因为样例不是随的，样例中原点正好在 \(V'\) 中（因为有球心直接落在了原点上！）。当然这个情况你可以自行随机一个向量把球心集体位移一下再位移回来。

#include <cmath>
#include <cstdio>
#include <algorithm>
using namespace std;
const int N=505;
int d,n;
inline double dot(double *f,double *g){
	double res=0;
	for(int i=1;i<=d;++i) res+=f[i]*g[i];
	return res;
}
double e[N][N],a[N][N],ra[N],rb[N],s[N][N];
int main(){
	scanf("%d%d",&d,&n);
	for(int i=1;i<=n;++i)
		for(int j=1;j<=d+1;++j) scanf("%lf",&e[i][j]);
	if(d<n){
		for(int i=1;i<=d+1;++i){
			for(int j=1;j<=d;++j) a[i][j]=e[i][j]*2;
			a[i][d+1]=-1;
			a[i][d+2]=dot(e[i],e[i])-e[i][d+1]*e[i][d+1];
		}
		for(int i=1;i<=d+1;++i){
			for(int j=d+2;j>=i;--j) a[i][j]/=a[i][i];
			for(int j=i+1;j<=d+1;++j)
				for(int k=d+2;k>=i;--k) a[j][k]-=a[j][i]*a[i][k];
		}
		for(int i=d+1;i;--i)
			for(int j=1;j<i;++j) a[j][d+2]-=a[j][i]*a[i][d+2];
		for(int i=1;i<=d;++i) printf("%.8lf ",a[i][d+2]);
		putchar('\n');
	}
	else{
		for(int i=1;i<=n;++i)
			for(int j=1;j<=i;++j) s[i][j]=dot(e[i],e[j]);
		for(int i=1;i<=n;++i)
			for(int j=i+1;j<=n;++j) s[i][j]=s[j][i];
		for(int i=1;i<=n;++i){
			for(int j=1;j<=n;++j) a[i][j]=s[i][j]*2;
			a[i][n+1]=s[i][i]-e[i][d+1]*e[i][d+1];
			a[i][n+2]=1;
		}
		for(int i=1;i<=n;++i){
			for(int j=n+2;j>=i;--j) a[i][j]/=a[i][i];
			for(int j=i+1;j<=n;++j)
				for(int k=n+2;k>=i;--k) a[j][k]-=a[j][i]*a[i][k];
		}
		for(int i=n;i;--i){
			for(int j=i-1;j;--j){
				a[j][n+1]-=a[j][i]*a[i][n+1];
				a[j][n+2]-=a[j][i]*a[i][n+2];
			}
			for(int j=1;j<=d;++j){
				ra[j]+=e[i][j]*a[i][n+1];
				rb[j]+=e[i][j]*a[i][n+2];
			}
		}
		double A=0,B=0,C=0;
		for(int i=1;i<=d;++i){
			A+=rb[i]*rb[i];
			B+=2*ra[i]*rb[i];
			C+=ra[i]*ra[i];
		}
		--B;
		double r=(-B+sqrt(B*B-4*A*C))/(2*A);
		for(int i=1;i<=d;++i) printf("%.8lf ",ra[i]+rb[i]*r);
		putchar('\n');
	}
	return 0;
}