P4035 [JSOI2008] 球形空间产生器

P4035 [JSOI2008] 球形空间产生器

大意

给你一个 \(n\) 维空间，给我 \(n + 1\) 个点，你要确定这个球的中心点的坐标。

思路

原始方程（\(n+1\) 个）对于每一个点 \(i\)，到球心 \((x_1, \dots, x_n)\) 的距离平方相等：$$\sum_{j=1}^n (p_{i,j} - x_j)^2 = R^2$$

降次消元（作差）用相邻两个方程 \((i+1)\) 和 \((i)\) 相减，消去二次项 \(x_j^2\) 和 \(R^2\)：$$\sum_{j=1}^n (p_{i+1,j} - x_j)^2 - \sum_{j=1}^n (p_{i,j} - x_j)^2 = 0$$展开并整理：$$\sum_{j=1}^n (p_{i+1,j}^2 - 2p_{i+1,j}x_j) - \sum_{j=1}^n (p_{i,j}^2 - 2p_{i,j}x_j) = 0$$

标准线性方程将未知数 \(x_j\) 移至左侧，常数移至右侧：$$\sum_{j=1}^n 2(p_{i+1,j} - p_{i,j})x_j = \sum_{j=1}^n (p_{i+1,j}^2 - p_{i,j}^2)$$

增广矩阵 \([A|B]\)

\[\begin{bmatrix} 2(p_{2,1}-p_{1,1}) & \dots & 2(p_{2,n}-p_{1,n}) & \mid & \sum (p_{2,j}^2 - p_{1,j}^2) \\ 2(p_{3,1}-p_{2,1}) & \dots & 2(p_{3,n}-p_{2,n}) & \mid & \sum (p_{3,j}^2 - p_{2,j}^2) \\ \vdots & \ddots & \vdots & \mid & \vdots \\ 2(p_{n+1,1}-p_{n,1}) & \dots & 2(p_{n+1,n}-p_{n,n}) & \mid & \sum (p_{n+1,j}^2 - p_{n,j}^2) \end{bmatrix} \]

代码

#include<iostream>
#include<cmath>
#include<iomanip>
using namespace std;

const int MAXN = 1005;
int n;
double xy[MAXN][MAXN];
double ans[MAXN];
double d[MAXN];
double nw[MAXN];

void Temper(){
    double tot = 0;
    for(int i = 1;i <= n + 1;i ++){
        d[i] = 0;
        nw[i] = 0;
        for(int j = 1;j <= n;j ++){
            d[i] += (xy[i][j] - ans[j]) * (xy[i][j] - ans[j]);
        }
        d[i] = sqrt(d[i]);
        tot += d[i];
    }
    tot /= (n + 1);
    for(int i = 1;i <= n + 1;i ++){
        for(int j = 1;j <= n;j ++){
            nw[j] += (d[i] - tot) * (xy[i][j] - ans[j]) / tot;
        }
    }
}

int main() {
    cin >> n;
    for(int i = 1;i <= n + 1;i ++){
        for(int j = 1;j <= n;j ++){
            cin >> xy[i][j];
            ans[j] += xy[i][j];
        }
    }
    for(int i = 1;i <= n;i ++){
        ans[i] /= (n + 1);
    }
    for(double t = 10001;t >= 0.0001;t *= 0.99995){
        Temper();
        for(int i = 1;i <= n;i ++){
            ans[i] += t * nw[i];
        }
    }
    cout << fixed << setprecision(3);
    for(int i = 1;i <= n;i ++){
        cout << ans[i] << ' ';
    }
    return 0;
}