洛谷 P4035：[JSOI2008] 球形空间产生器 ← 高斯-约旦消元法

​【题目来源】
https://www.luogu.com.cn/problem/P4035

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

【输入格式】
第一行是一个整数 n (1≤n≤10)。接下来的 n+1 行，每行有 n 个实数，表示球面上一点的 n 维坐标。每一个实数精确到小数点后 6 位，且其绝对值都不超过 20000。

【输出格式】
有且只有一行，依次给出球心的 n 维坐标（ n 个实数），两个实数之间用一个空格隔开。每个实数精确到小数点后 3 位。数据保证有解。你的答案必须和标准输出一模一样才能够得分。​​​​​​​

【输入样例】
2
0.0 0.0
-1.0 1.0
1.0 0.0​​​​​​​

【输出样例】
0.500 1.500

【数据范围】
1≤n≤10

【算法分析】
● 设球心 ，球面上两点 a,b，半径 R：
，
两式相减，得：
拿第 0 个点分别和第 1, 2, …, n 个点做差，得到一共 n 个方程，正好 n 元一次方程组，用高斯-约旦消元法求解。
本题保证唯一解，不需要判无解 / 无穷解。

● 高斯-约旦消元法应用
​​​​​​​洛谷 P3389：[模板] 高斯消元法：https://blog.csdn.net/hnjzsyjyj/article/details/162507345
AcWing 883：高斯消元解线性方程组：https://blog.csdn.net/hnjzsyjyj/article/details/170562311

【算法代码】

#include <bits/stdc++.h>
using namespace std;

const int N=15;
const double eps=1e-8;
double a[N][N];
double p[N][N];
int n;

void gauss() {
    int r=0;
    for(int c=0; c<n; c++) {
        //1.Select pivot
        int idx=r;
        for(int i=r; i<n; i++) {
            if(fabs(a[i][c])>fabs(a[idx][c])) idx=i;
        }
        if(fabs(a[idx][c])<eps) continue;
        swap(a[r],a[idx]);

        //2.Normalize the current row,
        //and the pivot becomes 1
        double div=a[r][c];
        for(int j=c; j<=n; j++) a[r][j]/=div;

        //3.Eliminate all rows
        for(int i=0; i<n; i++) {
            if(i!=r && fabs(a[i][c])>eps) {
                double val=a[i][c];
                for(int j=c; j<=n; j++) {
                    a[i][j]-=val*a[r][j];
                }
            }
        }
        r++;
    }
}

int main() {
    cin>>n;
    for(int i=0; i<=n; i++) {
        for(int j=0; j<n; j++) {
            cin>>p[i][j];
        }
    }

    for(int i=1; i<=n; i++) {
        double t=0;
        for(int j=0; j<n; j++) {
            a[i-1][j]=2*(p[0][j]-p[i][j]);
            t+=p[0][j]*p[0][j]-p[i][j]*p[i][j];
        }
        a[i-1][n]=t;
    }

    gauss();
    for(int i=0; i<n; i++) {
        printf("%.3lf ",a[i][n]);
    }
    return 0;
}

/*
in:
2
0.0 0.0
-1.0 1.0
1.0 0.0

out:
0.500 1.500
*/

【参考文献】
https://blog.csdn.net/hnjzsyjyj/article/details/162507345
https://blog.csdn.net/hnjzsyjyj/article/details/170562311

 

​