《算法竞赛进阶指南》0x35 高斯消元与线性空间

题目链接：https://www.acwing.com/problem/content/209/

给出n维空间中的n维圆上的n+1个点，要求求这个圆的圆心，保证一定有解，这样一共有n+1个n圆二次方程，差分之后变成n个n元方程，保证一定有解那么系数矩阵一定是满秩的，通过高斯消元可以求解

代码：

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<math.h>
using namespace std;
#define maxn 20
double a[maxn][maxn],b[maxn],c[maxn][maxn];
int n;
int main(){
    cin>>n;
    for(int i=1;i<=n+1;i++)
        for(int j=1;j<=n;j++)
            scanf("%lf",&a[i][j]);
    
    for(int i=1;i<=n;i++){    //构造n个差分方程，
        for(int j=1;j<=n;j++){
            c[i][j]=2*(a[i][j]-a[i+1][j]);
            b[i]+=a[i][j]*a[i][j]-a[i+1][j]*a[i+1][j];
        }
    }
    //高斯消元 
    for(int i=1;i<=n;i++){
        for(int j=i;j<=n;j++){//只交换比i大的方程j 
            if(fabs(c[j][i])>1e-8){//系数不为0 
                for(int k=1;k<=n;k++)swap(c[j][k],c[i][k]);//交换两个方程
                swap(b[i],b[j]); 
            }
        }
        
        for(int j=1;j<=n;j++){//第j行减去第i行的rate倍 
            if(j==i)continue;//将其他方程的对应项全部变成0 
            double rate=c[j][i]/c[i][i];
            for(int k=i;k<=n;k++){//i之前的全都是0，不用处理 
                c[j][k]-=c[i][k]*rate; 
            } 
            b[j]-=b[i]*rate;
        }
    }
    
    for(int i=1;i<=n;i++)b[i]/=c[i][i],c[i][i]=1;//将x_i的系数变为1， 
    for(int i=1;i<n;i++)printf("%.3lf ",b[i]);
    printf("%.3lf\n",b[n]);
}