[JSOI2008]球形空间产生器sphere

Description

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

Input

第一行是一个整数$n(1<=N=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

HINT

提示：
给出两个定义：
1、 球心：到球面上任意一点距离都相等的点。
2、 距离：设两个$n$为空间上的点$A, B$的坐标为$(a_1, a_2, …, a_n), (b_1, b_2, …, b_n)$，则$AB$的距离定义为：$dist = \sqrt{ (a_1-b_1)^2 + (a_2-b_2)^2 +… + (a_n-b_n)^2 }$

Source

思路

高斯消元模板。
设球心$O$的坐标为$(x_1,x_2,\cdots,x_n)$，点$P$在$j$方向上的坐标为$a_{P,j}$，则有：

$$\begin{align}\left|OA\right|&=(a_{A,1}-x_1)^2+(a_{A,2}-x_2)^2+(a_{A,3}-x_3)^2+\cdots+(a_{A,n}-x_n)^2\\&=a_{A,1}^2-2a_{A,1}x_1+x_1^2+\cdots+a_{A,n}^2+2a_{A,n}x_n^2+x_n^2\end{align}$$ 同理， $$\left|OB\right|=a_{B,1}^2-2a_{B,1}x_1+x_1^2+\cdots+a_{B,n}^2+2a_{B,n}x_n^2+x_n^2$$ 又因为 $$\left|OA\right|=\left|OB\right|$$ 所以 $$a_{A,1}^2-2a_{A,1}x_1+x_1^2+\cdots+a_{A,n}^2+2a_{A,n}x_n^2+x_n^2=a_{B,1}^2-2a_{B,1}x_1+x_1^2+\cdots+a_{B,n}^2+2a_{B,n}x_n^2+x_n^2$$ 整理得 $$\sum_{i=1}^n(2a_{A,i}-2a_{B,i})x_i=\sum_{i=1}^n(a_{A,i}^2-a_{B,i}^2)$$ 得出个方程，高斯消元解方程即可得出答案。

代码

#include <cstdio>
#include <cmath>
#include <algorithm>

const int maxn=10;
const double eps=0.0001;

int n;

struct matrix
{
  double a[maxn+1][maxn+2];

  inline int gauss()
  {
    for(register int  i=1; i<n; ++i)
      {
        int t=i;
        for(register int j=i+1; j<=n; ++j)
          {
            if(fabs(a[t][i])<fabs(a[j][i]))
              {
                t=j;
              }
          }
        if(fabs(a[t][i])<eps)
          {
            return 1;
          }
        for(register int j=i; j<=n+1; ++j)
          {
            std::swap(a[t][j],a[i][j]);
          }
        for(register int j=i+1; j<=n; ++j)
          {
            if(!a[j][i])
              {
                continue;
              }
            double f=a[j][i]/a[i][i];
            for(register int k=i; k<=n+1; ++k)
              {
                a[j][k]-=f*a[i][k];
              }
          }
      }
    for(register int i=n; i; --i)
      {
        a[i][n+1]=a[i][n+1]/a[i][i];
        a[i][i]=1;
        for(register int j=1; j<i; ++j)
          {
            a[j][n+1]-=a[j][i]*a[i][n+1];
          }
      }
    return 0;
  }

  inline int print_ans()
  {
    for(register int i=1; i<n; ++i)
      {
        printf("%.3lf ",a[i][n+1]);
      }
    printf("%.3lf",a[n][n+1]);
    return 0;
  }
};

matrix t;
double point[maxn+2][maxn+1];

int main()
{
  scanf("%d",&n);
  for(register int i=1; i<=n+1; ++i)
    {
      for(register int j=1; j<=n; ++j)
        {
          scanf("%lf",&point[i][j]);
        }
    }
  for(register int i=1; i<=n; ++i)
    {
      double sum=0;
      for(register int j=1; j<=n; ++j)
        {
          t.a[i][j]=2*(point[i][j]-point[i+1][j]);
          sum+=point[i][j]*point[i][j]-point[i+1][j]*point[i+1][j];
        }
      t.a[i][n+1]=sum;
    }
  t.gauss();
  t.print_ans();
  return 0;
}