球形空间产生器 bzoj 1013 高斯消元

题目大意

给出n维空间中给出n+1个点的坐标，求出球心坐标。

分析

设球心坐标为（x1,x2,x3…xn）
那么就有
(a1-x1)^2+(a2-x2)^2+…(an-xn)^2=r^2
(b1-x1)^2+(b2-x2)^2+…(bn-xn)^2=r^2
….
只要拿后n个方程分别去减第一个方程，就可以得到n个一次方程了
2*(a1-b1)x1+2*(a2-b2)x2+…..+2(an-bn)xn=a1^2-b1^2+a2^2-b2^2….an^2-bn^2
…
高斯消元即可

code

#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>

using namespace std;

double x[1000];
int n;

double b[1000];
double a[1000][1000];

double y[1000];
int main()
{
    scanf("%d",&n);
    for (int j=1;j<=n;j++)
        cin>>y[j];
    for (int i=2;i<=n+1;i++)
    {
        for (int j=1;j<=n;j++)
        {
            double xx;
            cin>>xx;
            a[i-1][j]=2*(xx-y[j]);
            b[i-1]+=(xx*xx)-(y[j]*y[j]);
        }
    }   
    for (int i=1;i<=n;i++)
    {
        double maxx=fabs(a[i][i]);
        int k=i;
        double c=0;
        for (int j=i+1;j<=n;j++)
            if (fabs(a[j][i])>maxx)
            {
                maxx=fabs(a[j][i]);
                k=j;
            }
        for (int j=1;j<=n;j++)
        {
            c=a[i][j];
            a[i][j]=a[k][j];
            a[k][j]=c;
        }
        c=b[i]; b[i]=b[k]; b[k]=c;
        for (int j=1;j<=n;j++)
        {
            if ((j==i)||(!a[j][i])) continue;\
            double t=a[i][i]/a[j][i];
            for (int l=1;l<=n;l++)
                a[i][l]/=t;
            b[i]/=t;
            b[j]-=b[i];
            for (int l=1;l<=n;l++)
                a[j][l]-=a[i][l];
        }
    }
    for (int i=1;i<=n;i++)
        x[i]=b[i]/a[i][i];
    for (int i=1;i<n;i++)
    {
        printf("%.3f ",x[i]);
    }
    printf("%.3f",x[n]);   
}