HYSBZ 1013 球形空间产生器sphere（高斯消元）

给定n维空间下的n+1个点，求这n个点所在的球面的球心

我们设球心为X(x1,x2,...,xn)

假设有两点A(a1,a2,...,an)和B(b1,b2,...,bn)

那么我们可以得到两个方程

(x1-a1)^2+(x2-a2)^2+...+(xn-an)^2=r^2

(x1-b1)^2+(x2-b2)^2+...+(xn-bn)^2=r^2

这些方程都是二次的，无法套用高斯消元

但是我们可以做一些处理 将上面两个方程相减可得

(a1-b1)x1+(a2-b2)x2+...+(an-bn)xn=[ (a1^2-b1^2)+(a2^2-b2^2)+...+(an^2-bn^2) ]/2

r被消掉，n个方程，n个未知数套用高斯消元模板即可

# include <stdio.h>
# include <algorithm>
# include <string.h>
# include <math.h>
using namespace std;
#define eps 1e-9
const int MAXN=220;
double a[MAXN][MAXN],x[MAXN];//方程的左边的矩阵和等式右边的值，求解之后x存的就是结果
int equ,var;//方程数和未知数个数
/*
*返回0表示无解，1表示有解
*/
void Gauss()
{
    int i,j,k,col,max_r;
    for(k=0,col=0; k<equ&&col<var; k++,col++)
    {
        max_r=k;
        for(i=k+1; i<equ; i++)
            if(fabs(a[i][col])>fabs(a[max_r][col]))
                max_r=i;
        if(fabs(a[max_r][col])<eps)
            return ;
        if(k!=max_r)
        {
            for(j=col; j<var; j++)
                swap(a[k][j],a[max_r][j]);
            swap(x[k],x[max_r]);
        }
        x[k]/=a[k][col];
        for(j=col+1; j<var; j++)a[k][j]/=a[k][col];
        a[k][col]=1;
        for(i=0; i<equ; i++)
            if(i!=k)
            {
                x[i]-=x[k]*a[i][k];
                for(j=col+1; j<var; j++)a[i][j]-=a[k][j]*a[i][col];
                a[i][col]=0;
            }
    }
    return ;

}
///(a1-b1)x1+(a2-b2)x2+....+(an-bn)x = [(a1^2-b1^2)+(a2^2-b2^2)+...+(an^2-bn^2)]/2

int main()
{
    int n;
    double p[15][15];
    while(~scanf("%d",&n))
    {
        equ=var=n;
        for(int i = 0; i <= n; i++)
            for(int j = 0; j < n; j++)
                scanf("%lf", &p[i][j]);
        for(int i = 0; i < n; i++)
            {
                for(int j = 0; j < n; j++)
            {
                a[i][j] = (p[i][j] - p[i+1][j]);
                x[i] += p[i][j] * p[i][j];
                x[i] -= p[i+1][j] * p[i+1][j];
            }
            x[i]/=2;
            }
           Gauss();
           for(int i=0;i<n;i++)
           {
               if(i==n-1)
                printf("%.3lf\n",x[i]);
               else
                printf("%.3lf ",x[i]);
           }

    }
    return 0;
}