洛谷 P4035 [JSOI2008]球形空间产生器
思路
高斯消元
题意：在(n)维的球形空间中给定(n+1)个点，求到所有(n+1)个点的距离相等的点的坐标
由题意易知我们要求出在(n)维空间中的一个点((x_1,x_2,x_3,…x_n))，满足：
[forall iin[1,n+ 1],sumlimits_{j = 1}^{n}(a_{i,j} - x_j)^2 = R
]
其中(R)是一个常数，第(i)个点的坐标为((a_{i,1},a_{i,2},a_{i,3},…a_{i,n}))
假设有(i_1,i2in[1,n+1],i1≠i2)，由(sumlimits_{j = 1}^{n}(a_{i_1,j} - x_j)^2 = R)和(sumlimits_{j = 1}^{n}(a_{i_2,j} - x_j)^2 = R)得
[sumlimits_{j = 1}^{n}(a_{i_1,j} - x_j)^2 = sumlimits_{j = 1}^{n}(a_{i_2,j} - x_j)^2
]
展开式子得
[sumlimits_{j = 1}^{n}(a_{i_1,j}^2 + x_j^2 - 2 *(a_{i_1,j}* x_j))=sumlimits_{j = 1}^{n}(a_{i_2,j}^2 + x_j^2 - 2 *(a_{i_2,j}* x_j))
]
进一步化简
[sumlimits_{j = 1}^{n}(a_{i_1,j}^2 - 2 *(a_{i_1,j}* x_j))=sumlimits_{j = 1}^{n}(a_{i_2,j}^2 - 2 *(a_{i_2,j}* x_j))
]
[sumlimits_{j = 1}^{n} 2 *(a_{i_1,j} - a_{i_2,j})x_j= sumlimits_{j = 1}^{n}(a_{i_1,j}^2 - a_{i_2,j}^2)
]
这样就转化成了一个线性方程组，由此可以用每个(i)和(i+1)两两组合，得出(n)个线性方程组，则最后的矩阵为
[egin{bmatrix}2(a_{1,1}-a_{2,1}) 2(a_{1,2}-a_{2,2}) … 2(a_{1,n}-a_{2,n})sumlimits_{j=1}^{n}(a^2_{1,j}-a^2_{2,j})\2(a_{2,1}-a_{3,1}) 2(a_{2,2}-a_{3,2}) … 2(a_{2,n}-a_{3,n})sumlimits_{j=1}^{n}(a^2_{2,j}-a^2_{3,j})\2(a_{3,1}-a_{4,1}) 2(a_{3,2}-a_{4,2}) … 2(a_{3,n}-a_{4,n})sumlimits_{j=1}^{n}(a^2_{3,j}-a^2_{4,j})\…\2(a_{n,1}-a_{n+1,1}) 2(a_{n,2}-a_{n+1,2}) … 2(a_{n,n}-a_{n+1,n})sumlimits_{j=1}^{n}(a^2_{n,j}-a^2_{n+1,j})end{bmatrix}
]
对此矩阵进行高斯消元求解即可，由于保证有解，所以直接做就好了
代码
/*
Author:loceaner
*/
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define eps 1e-8
using namespace std;
const int A = 22;
const int B = 1e6 + 11;
const int mod = 1e9 + 7;
const int inf = 0x3f3f3f3f;
inline int read() {
char c = getchar(); int x = 0, f = 1;
for( ; !isdigit(c); c = getchar()) if(c == '-') f = -1;
for( ; isdigit(c); c = getchar()) x = x * 10 + (c ^ 48);
return x * f;
}
int n;
double a[A][A], G[A][A], b[A];
//a是输入的点坐标，G是矩阵序数数组，b是方程右边的常数
int main() {
n = read();
for (int i = 1; i <= n + 1; i++)
for (int j = 1; j <= n; j++) scanf("%lf", &a[i][j]);
for (int i = 1; i <= n; i++)
for (int j = 1; j <= n; j++) {
G[i][j] = 2 * (a[i][j] - a[i + 1][j]);
b[i] = b[i] + a[i][j] * a[i][j] - a[i + 1][j] * a[i + 1][j];
}
for (int i = 1; i <= n; i++) {
for (int j = i; j <= n; j++) {
if (fabs(G[i][j]) > eps) {
for (int k = 1; k <= n; k++) swap(G[i][k], G[j][k]);
swap(b[i], b[j]);
}
}
for (int j = 1; j <= n; j++) {
if (i == j) continue;
double tmp = G[j][i] / G[i][i];
for (int k = i; k <= n; k++) G[j][k] -= G[i][k] * tmp;
b[j] -= b[i] * tmp;
}
}
for (int i = 1; i <= n; i++) printf("%.3lf ", b[i] / G[i][i]);
return 0;
}