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

$球型空间产生器$球型空间产生器

$Description$Description

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

$Solution$Solution

题意: 求 $n$n维点 $O$O( $a_1,a_2,...,a_n$a1​,a2​,...,an​), 使 $O$O 距离 球面上 $n+1$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$ (a1​−b1,1​)2+(a2​−b1,2​)2+...+(an​−b1,n​)2=D2 (a1​−b2,1​)2+(a2​−b2,2​)2+...+(an​−b2,n​)2=D2... (a1​−bn+1,1​)2+(a2​−bn+1,2​)2+...+(an​−bn+1,n​)2=D2 

相邻式子相减消去 $D^2$D2, 得到
$(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..略$ (2a1​−b1,1​−b2,1​)(b2,1​−b1,1​)+...+(2an​−b1,n​−b2,n​)(b2,n​−b1,n​)=0..略 

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

系数 $A[i,j]=2\ast \left(b\right[i+1,j]-b[i,j\left]\right)$A[i,j]=2∗(b[i+1,j]−b[i,j])

常数 $A[i,n+1]=\sum _{j=1}^n\left(b\right[i,j]+b[i+1,j\left]\right)\ast \left(b\right[i+1,j]-b[i,j\left]\right)$A[i,n+1]=j=1∑n​ (b[i,j]+b[i+1,j]) ∗ (b[i+1,j]−b[i,j])

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

$Code$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;
}