[JSOI2008]球形空间产生器（线性代数+高斯消元）

Posted 2020-12-11 _CRしち_のblog

题目大意

给你一个n维球体上的n+1个点，让你求这个n维球体的球心。数据保证球心是唯一的。

Analysis

将球心设出来为$(x_1, x_2, cdots, x_n)$，设半径为$r$。设球上一点为$(y_1, y_2, cdots, y_n)$，根据n维空间内两点之间距离公式得$sum_{i=1}^n (y_i-x_i)^2=r^2$。

设每个点表示为$(a_{i,1}, a_{i,2}, cdots, a_{i,n})$。可以得到方程组：

$$egin{cases}
sum_{i=1}^n (a_{1,i}-x_i)^2=r^2\
sum_{i=1}^n (a_{2,i}-x_i)^2=r^2\
cdots\
sum_{i=1}^n (a_{n+1,i}-x_i)^2=r^2\
end{cases}
$$

将第$i$个方程与第$i+1$个方程相减得：

$$egin{cases}
sum_{i=1}^n (a_{1,i}^2-a_{2,i}^2-2x_i(a_{1,i}-a_{2,i}))=0\
sum_{i=1}^n (a_{2,i}^2-a_{3,i}^2-2x_i(a_{2,i}-a_{3,i}))=0\
cdots\
sum_{i=1}^n (a_{n,i}^2-a_{n+1,i}^2-2x_i(a_{n,i}-a_{n+1,i}))=0\
end{cases}
$$

再变成标准形式：

$$egin{cases}
sum_{i=1}^n (2x_i(a_{1,i}-a_{2,i}))=sum_{i=1}^n a_{1,i}^2-a_{2,i}^2\
sum_{i=1}^n (2x_i(a_{2,i}-a_{3,i}))=sum_{i=1}^n a_{2,i}^2-a_{3,i}^2\
cdots\
sum_{i=1}^n (2x_i(a_{n,i}-a_{n+1,i}))=sum_{i=1}^n a_{n,i}^2-a_{n+1,i}^2\
end{cases}
$$

n个未知数n个方程，直接高斯消元。因为题目保证有解所以不需要做过多的判断。

Code：

#include <iostream>
#include <cstdio>
using namespace std;
const int N = 11;
int n;
double a[N][N], A[N][N], B[N];
void Gauss() {
    
    for (int i = 1; i <= n; i++) {
        int mx = i;
        for (int j = i + 1; j <= n; j++) {
            if (A[j][i] > A[mx][i]) {
                mx = j;
            }
        }
        for (int j = i; j <= n; j++) {
            swap(A[i][j], A[mx][j]);
        }
        swap(B[i], B[mx]);
        if (!A[i][i]) continue;
        for (int j = 1; j <= n; j++) {
            if (i == j) continue;
            double Rate = A[j][i] / A[i][i];
            for (int k = i; k <= n; k++) {
                A[j][k] -= A[i][k] * Rate;
            }
            B[j] -= B[i] * Rate;
        }
    }
}
int main() {
    scanf("%d", &n);
    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++) {
            A[i][j] = 2.0 * (a[i][j] - a[i + 1][j]);
            B[i] += a[i][j] * a[i][j] - a[i + 1][j] * a[i + 1][j];
        }
    }
    Gauss();
    for (int i = 1; i <= n; i++) printf("%.3lf ", B[i] / A[i][i]);
    return 0;
}