球形空间产生器 [高斯消元/模拟退火]

$球型空间产生器$

---

$\mathcal{Description}$

有一个球形空间产生器能够在 $n$ 维空间中产生一个坚硬的球体。现在，你被困在了这个 n 维球体中，你只知道球面上 $n+1$ 个点的坐标，你需要以最快的速度确定这个 $n$ 维球体的球心坐标，以便于摧毁这个球形空间产生器。
$n <= 10$

---

$\mathcal{Solution}$

题意: 求 $n$维点 $O$($a_1,a_2,...,a_n$), 使 $O$ 距离 球面上 $n+1$ 个点的 距离相等,

使用 式子 表示出来, 即
$\ \\ (a_1-b_{1,1})^2+(a_2-b_{1,2})^2+...+(a_n-b_{1,n})^2=D^2\\ \ \\ (a_1-b_{2,1})^2+(a_2-b_{2,2})^2+...+(a_n-b_{2,n})^2=D^2\\ . \\ . \\ . \\ \ \\ (a_1-b_{n+1,1})^2+(a_2-b_{n+1,2})^2+...+(a_n-b_{n+1,n})^2=D^2 \\ \ \\$

相邻式子相减消去 $D^2$, 得到
$\ \\ (2a_1-b_{1,1}-b_{2,1})(b_{2,1}-b_{1,1})+...+(2a_n-b_{1,n}-b_{2,n})(b_{2,n}-b_{1,n})=0\\ . \\ . \\ 略 \\ \ \\$

发现已经是一个 线性方程组 了,

系数 $A[i,j]=2*(b[i+1,j]-b[i,j])$

常数 $A[i,n+1]=\sum_{j=1}^n \ (b[i,j] + b[i+1,j])\ *\ (b[i+1,j]-b[i,j])$

于是 高斯消元 $O(N^3)$ 求解 .

---

$\mathcal{Code}$

#include<bits/stdc++.h>
#define reg register

const int maxn = 15;

int N;
double A[maxn][maxn];
double B[maxn][maxn];

int main(){
        scanf("%d", &N);
        for(reg int i = 1; i <= N+1; i ++)
                for(reg int j = 1; j <= N; j ++) scanf("%lf", &B[i][j]);
        for(reg int i = 1; i <= N; i ++){
                for(reg int j = 1; j <= N; j ++) A[i][j] = 2.0*(B[i+1][j] - B[i][j]);
                A[i][N+1] = 0;
                for(reg int j = 1; j <= N; j ++)
                        A[i][N+1] += (B[i][j]+B[i+1][j])*(B[i+1][j] - B[i][j]);
        } 
        for(reg int i = 1; i <= N; i ++){
                int max_id = i;
                for(reg int j = i+1; j <= N; j ++)
                        if(fabs(A[max_id][i]) < fabs(A[j][i])) max_id = j;
                std::swap(A[i], A[max_id]);
                double tmp = A[i][i];
                for(reg int j = i; j <= N+1; j ++) A[i][j] /= tmp;
                for(reg int j = i+1; j <= N; j ++){
                        tmp = A[j][i];
                        for(reg int k = i; k <= N+1; k ++)
                                A[j][k] -= A[i][k] * tmp;
                }
        }
        for(reg int i = N; i >= 1; i --)
                for(reg int j = i+1; j <= N; j ++) A[i][N+1] -= A[i][j]*A[j][N+1];
        for(reg int i = 1; i <= N; i ++) printf("%.3lf ", A[i][N+1]);
        return 0;
}