算法竞赛——进阶指南——acwing207. 球形空间产生器 高斯消元应用

n维球体的通式：

(x1-a1)^2 + (x2-a2)^2 + …… + (xn-an)^2  =R ^ 2

球心坐标就是(a1,a2,……,an)

给出n+1个坐标，均带入，发现不好解，因为是二次方程，而且R^2是未知数

由于所有方程右边都是R^2  我们可以两两联立方程，消去一次 变成多元一次方程 用高斯消元搞搞就行。

例如：

(x11-a1)^2 + (x12-a2)^2 + …… + (x12-an)^2  =R ^ 2

(x21-a1)^2 + (x22-a2)^2 + …… + (x2n-an)^2  =R ^ 2

联立前两个得到：

2*(x21-x11)*a1+2*(x22-x12)+……+ 2 * (x2n-x1n) =(x21 ^ 2 - x11 ^ 2)+(x22 ^ 2 - x12 ^ 2)+……+(x2n ^ 2 - x1n ^ 2)

刚好是高斯消元初始矩阵的一行

两两联立即可

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define ls (o<<1)
#define rs (o<<1|1)
#define pb push_back
const double PI= acos(-1.0);
const double eps=1e-8;
const int M = 20;
double a[M][M],b[M][M];
int n;
int main()
{
    scanf("%d",&n);
    for(int i=1;i<=n+1;i++)
    for(int j=1;j<=n;j++)scanf("%lf",&b[i][j]);
    
    for(int i=1;i<=n;++i)
        for(int j=1;j<=n;++j)
        	a[i][j]=2*(b[i+1][j]-b[i][j]),a[i][n+1]+=b[i+1][j]*b[i+1][j]-b[i][j]*b[i][j];
    //下面高斯消元解方程 
    for(int i=1;i<=n;++i)//枚举列（项） 
    {
        int max=i;
        for(int j=i+1;j<=n;++j)//选出该列最大系数 
            if(fabs(a[j][i])>fabs(a[max][i]))max=j;
    	if(i!=max)swap(a[i],a[max]);//交换
        if(fabs(a[i][i])<eps)//最大值等于0则说明该列都为0，肯定无解 
        {
            puts("No Solution");
            return 0;
        }
        for(int j=1;j<=n;++j)//每一项都减去一个数（就是小学加减消元）
        {
            if(j==i)continue;
            double temp=a[j][i]/a[i][i];
            for(int k=i+1;k<=n+1;++k)
                a[j][k]-=a[i][k]*temp;
        }
    }
    for(int i=1;i<=n;++i)
        printf("%.3lf ",a[i][n+1]/a[i][i]);
    return 0;
}