本文链接：https://blog.csdn.net/weixin_43960414/article/details/120856728

题面

有一个 $n$ 维球，给定 $n + 1$ 个在球面上的点，求球心坐标。

$n\leq 10$ 。

题解

好久以前的题了，昨天首 A 。

$n$ 太小了！明明可以开 100！

看到题解里有两种做法：

模拟退火
距离半径法高斯消元

后面那种，以前一直搞不懂，不知道怎么解二次的。于是很久都没做出来。

现在，我有一种更好的方法：向量高斯消元法。

我们知道，对于球上两个点 $A, B$ 和球心 $P$ ，满足： $P$ 到 $A B$ 中点的连线与 $A B$ 垂直。

于是，我们就可以列出这样的等式：
$\stackrel{\longrightarrow}{AB}\cdot(\stackrel{\longrightarrow}{AP}+\stackrel{\longrightarrow}{BP})=0\\$

即
$(\stackrel{\longrightarrow}{OB}-\stackrel{\longrightarrow}{OA})\cdot(2\stackrel{\longrightarrow}{OP}-(\stackrel{\longrightarrow}{OB}+\stackrel{\longrightarrow}{OA}))=0\\ (\stackrel{\longrightarrow}{OB}-\stackrel{\longrightarrow}{OA})\cdot2\stackrel{\longrightarrow}{OP}=(\stackrel{\longrightarrow}{OB}-\stackrel{\longrightarrow}{OA})\cdot(\stackrel{\longrightarrow}{OB}+\stackrel{\longrightarrow}{OA})=\stackrel{\longrightarrow}{OB}^2-\stackrel{\longrightarrow}{OA}^2$

从坐标来看，就是
$\sum_{i=1}^n2(X_B[i]-X_A[i])\cdot \underline{~X_P[i]~}=|OB|^2-|OA|^2$

其中画横线的就是方程未知量，右边是两个距离的平方相减。

我们把 $n + 1$ 个点中每两个相邻的点列一个这样的方程。

下一步拉板子。

CODE

一定有解，所以不考虑这么多了。

#include<set>
#include<map>
#include<queue>
#include<stack>
#include<vector>
#include<bitset>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
#define MAXN 205
#define LL long long
#define DB double
#define lowbit(x) (-(x) & (x))
#define ENDL putchar('\n')
#define FI first
#define SE second
#define SQ 450
#define eps 1e-9
LL read() {
    LL f=1,x=0;int s = getchar(); if(s<0) return -1;
    while(s < '0' || s > '9') {if(s=='-')f = -f;s = getchar();}
    while(s >= '0' && s <= '9') {x=x*10+(s^48);s = getchar();}
    return f*x;
}
void putpos(LL x) {if(!x)return ;putpos(x/10);putchar('0'+(x%10));}
void putnum(LL x) {
    if(!x) {putchar('0');return ;}
    if(x<0) {putchar('-');x = -x;}
    return putpos(x);
}
void AIput(LL x,int c) {putnum(x);putchar(c);}

int MOD = 1;
int n,m,s,o,k;
DB pt[105][105];
DB a[105][105];
DB Abs(DB x) {return x<0 ? -x:x;}

void Gauss(int n) {
	for(int i = 1;i <= n;i ++) {
		for(int j = i+1;j <= n;j ++) {
			if(Abs(a[j][i]) > Abs(a[i][i])) {
				swap(a[j],a[i]); break;
			}
		}
		for(int j = i+1;j <= n;j ++) {
			DB nm = a[j][i] / a[i][i];
			for(int k = n+1;k >= i;k --) {
				a[j][k] -= a[i][k] * nm;
			}
		}
	}
	for(int i = n;i > 0;i --) {
		for(int j = n;j > i;j --) {
			a[i][n+1] -= a[i][j] * a[j][n+1];
		}
		a[i][n+1] /= a[i][i];
	}return ;
}
int main() {
	n = read();
	for(int i = 1;i <= n+1;i ++) {
		for(int j = 1;j <= n;j ++) {
			scanf("%lf",&pt[i][j]);
		}
	}
	for(int i = 1;i <= n;i ++) {
		for(int j = 1;j <= n;j ++) {
			a[i][j] = 2.0*(pt[i+1][j]-pt[i][j]);
			a[i][n+1] += pt[i+1][j]*pt[i+1][j] - pt[i][j]*pt[i][j];
		}
	}
	Gauss(n);
	for(int i = 1;i <= n;i ++) {
		printf("%.3f",a[i][n+1]);
		putchar(i==n ? '\n':' ');
	}
    return 0;
}