题意
给出(n)维直角坐标系中(n + 1)个点的坐标，它们都在一个(n)维球面上，求球心坐标。
题解
设球面上某两个点坐标为((a_1, a_2, ... a_n))和((b_1, b_2, ... b_n))，则可以列出方程:
[(x_1 - a_1)^2 + (x_2 - a_2)^2 + ... + (x_n - a_n)^2 = (x_1 - b_1)^2 + (x_2 - b_2)^2 + ... + (x_n - b_n)^2
]
括号打开化简得
[2*(a_1 - b_1)x_1 + 2*(a_2 - b_2)x_2 + ... + 2*(a_n - b_n)x_n = a_1^2 - b_1^2 + a_2^2 - b_2^2 + ... + a_n^2 - b_n^2
]
这样可得n个方程，然后高斯消元即可。
#include <cstdio>
#include <cmath>
#include <cstring>
#include <algorithm>
#define space putchar(' ')
#define enter putchar('
')
using namespace std;
typedef long long ll;
template <class T>
void read(T &x){
char c;
bool op = 0;
while(c = getchar(), c < '0' || c > '9')
if(c == '-') op = 1;
x = c - '0';
while(c = getchar(), c >= '0' && c <= '9')
x = x * 10 + c - '0';
if(op) x = -x;
}
template <class T>
void write(T x){
if(x < 0) putchar('-'), x = -x;
if(x >= 10) write(x / 10);
putchar('0' + x % 10);
}
const int N = 20;
int n;
double g[N][N], f[N][N];
int main(){
/*
下面以样例输入为例，模拟一下高斯消元的过程。
样例输入：
2
0.0 0.0
-1.0 1.0
1.0 0.0
*/
scanf("%d", &n);
for(int i = 1; i <= n + 1; i++)
for(int j = 1; j <= n; j++)
scanf("%lf", &g[i][j]);
//首先，我们构造一个矩阵来表示n个方程。
for(int i = 1; i <= n; i++)
for(int j = 1; j <= n; j++){
f[i][j] = 2 * (g[i][j] - g[i + 1][j]);
f[i][n + 1] += g[i][j] * g[i][j] - g[i + 1][j] * g[i + 1][j];
}
/*
上一步中，我们构造出了n个方程，分别是：
2 * x[1] - 2 * x[2] = -2
-4 * x[1] + 2 * x[2] = 1
写成矩阵的形式就是：
2 -2 -2
-4 2 1
这个矩阵存在f[][]数组中。
*/
for(int i = 1; i <= n; i++){
//这次循环，我们以x[i]为主元。
int l = i;
for(int j = i + 1; j <= n; j++)
if(fabs(f[j][i]) > fabs(f[l][i])) l = j;
if(l != i)
for(int j = i; j <= n + 1; j++)
swap(f[i][j], f[l][j]);
//上一步中，我们找出了主元系数绝对值最大的一个方程，把它换到第i行（据说这么做精度能高一些）
/*
程序第一次执行到这里的时候，矩阵变成了
-4 2 1
2 -2 -2
*/
for(int j = n + 1; j >= i; j--) //因为循环内要用到当前的f[i][i]，所以f[i][i]最后修改
f[i][j] /= f[i][i]; //将主元系数变为1，方程中其他项系数也等比例扩大/缩小。
for(int j = i + 1; j <= n; j++)
for(int k = n + 1; k >= i; k--) //因为循环内要用到当前的f[j][i]，所以f[j][i]最后修改
f[j][k] -= f[i][k] * f[j][i];
/*
将第i个方程的每项系数扩大/缩小，使得主元系数和第j个方程主元系数相同
然后第j个方程 -= 第i个方程，这样第j个方程就消去了当前主元。
*/
/*
程序第一次进行到这一步的时候，矩阵变成了
1 -0.5 -0.25
0 -1 -1.5
第二次则变成了
1 -0.5 -0.25
0 1 1.5
*/
}
for(int i = n; i; i--) //循环到i的时候，f[i][n + 1]表示的已经是最后的解x[i]
for(int j = 1; j < i; j++)
f[j][n + 1] -= f[j][i] * f[i][n + 1]; //将x[i]带入到第j个方程中
for(int i = 1; i <= n; i++)
printf("%.3lf%c", f[i][n + 1], "
"[i == n]);
return 0;
}