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

Description

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

Input

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

Output

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

Sample Input

2

0.0 0.0

-1.0 1.0

1.0 0.0

Sample Output

0.500 1.500

Solution

记第$i$个点的坐标为$(x_{i,1},x_{i,2},...,x_{i,n})$,$0\le i\le n$，球心坐标为$(x_1,x_2,...,x_n)$

则由这$n+1$个点到球心距离相同知$\sum\limits_{j=1}^n (x_{i,j}-x_j)^2=\sum\limits_{j=1}^n (x_{0,j}-x_j)^2,i=1,2,...,n$

展开有$\sum\limits_{j=1}^n 2(x_{i,j}-x_{0,j})x_j=\sum\limits_{j=1}^n (x_{i,j}^2-x_{0,j}^2)$

用高斯消元解出$x_1,...,x_n$即可

Code

#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<vector>
#include<queue>
#include<map>
#include<set>
#include<ctime>
using namespace std;
typedef long long ll;
typedef pair<int,int>P;
const int INF=0x3f3f3f3f,maxn=15;
#define eps 1e-8
int sign(double x)
{
    if(fabs(x)<eps)return 0;
    if(x>eps)return 1;
    return -1;
}
double A[maxn][maxn],x[maxn];
int Gauss(int equ,int var)
{
    int i=1;
    for(;i<=equ;i++)
    {
        int pos=i;
        for(int j=i+1;j<=equ;j++)
            if(sign(abs(A[j][i])-abs(A[pos][i]))>0)pos=j;
        if(pos!=i)
            for(int j=i;j<=var+1;j++)swap(A[i][j],A[pos][j]);
        if(sign(A[i][i])==0)
        {
            i--;
            continue;
        }
        for(int j=i+1;j<=var+1;j++)A[i][j]/=A[i][i];
        A[i][i]=1;
        for(int j=i+1;j<=equ;j++)
            if(sign(A[j][i]))
            {
                for(int k=i+1;k<=var+1;k++)A[j][k]-=A[j][i]*A[i][k];
                A[j][i]=0;
            }
    }
    for(int j=i;j<=equ;j++)
        if(sign(A[j][var+1]))return -1;//无解
    if(i<=var)return 0;//无穷解，自由变元var-i+1个
    for(int i=var;i>=1;i--)
    {
        double temp=A[i][var+1];
        for(int j=i+1;j<=var;j++)temp-=A[i][j]*x[j];
        x[i]=temp/A[i][i];
    }
    return 1;
}
int n;
double a[maxn];
int main()
{
    scanf("%d",&n);
    for(int i=1;i<=n;i++)scanf("%lf",&a[i]);
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
        {
            double temp;
            scanf("%lf",&temp);
            A[i][j]=2.0*(temp-a[j]);
            A[i][n+1]+=temp*temp-a[j]*a[j];
        }
    Gauss(n,n);
    for(int i=1;i<=n;i++)printf("%.3f%c",x[i],i==n?'\n':' ');
    return 0;
}