BZOJ 1013 [JSOI2008]球形空间产生器sphere 【高斯消元】

Description

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

HINT

1<=n<=10

提示：给出两个定义：
1、 球心：到球面上任意一点距离都相等的点。
2、 距离：设两个n为空间上的点A, B的坐标为(a1, a2, …, an), (b1, b2, …, bn)，则AB的距离定义为：dist = sqrt( (a1-b1)^2 + (a2-b2)^2 + … + (an-bn)^2 )

Solution

题意：解 n 个 n 元一次方程

题目比较简单，主要是学一下高斯消元

根据定义很容易写出 n+1 个方程

发现对于每一个方程都有 n 个系数相等的二次项，

于是每两个方程做一下差

得到新的n个一次方程

当 n 等于 2 的时候大约是长这样的

2(a1-a2)x+2(b1-b2)y=a1^2-a2^2+b1^2-b2^2

矩阵中只保留各项系数

得到一个 n * ( n+1 ) 的矩阵，n+1 项表示其方程的值

对于矩阵上 f[ i ][ i ] 一项，我们尝试消去除 i 之后的方程上的第 i 位的系数，也就是方程做差消掉一个未知数，保留第 i 个方程对应第 i 个未知数

对于每一个 i  处理完之后，此时得到的矩阵所有的值分布在左上右下对角线的上方

然后依次拿下面的已知项将上方的方程中对应项的系数消掉

如果不考虑多解或无解的情况，最后可以解出一个对角线为1其他部分为0的矩阵，并且在 n+1 位置上上有各各未知数的唯一解

可以看这个例子感受一下：

 

2 1 -1 8

-3 -1 2 -11

-2 1 2 -3

 

跟着以上的方法来运算，这个矩阵可以转变为以下的样子：

 

2 1 -1 8

0 1/2 1/2 1

0 0 -1 1

 

这矩阵叫做“行梯列式”。

最后，可以利用同样的算法产生以下的矩阵，便可把所得出的解或其限制简明地表示出来：

 

1 0 0 2

0 1 0 3

0 0 1 -1

 

最后这矩阵叫做“简化行梯列式”，亦是高斯－约当消元法指定的步骤

复杂度O ( n^3 )

由于题目保证有解，随便写写就过了，然而喜闻乐见地PE了

#include<bits/stdc++.h>

#define maxn 10+5
#define set(a,b) memset(a,(b),sizeof(a))
#define fr(i,a,b) for(ll i=(a),_end_=(b);i<=_end_;i++)
#define rf(i,b,a) for(ll i=(a),_end_=(b);i>=_end_;i--)

using namespace std;

typedef long long ll;

double f[maxn][maxn],st[maxn][maxn];
int n,m;

double diff(double x,double y)
{
  return x*x-y*y;
}

void print()
{
  fr(i,1,n){
    fr(j,1,n+1)
      cout << f[i][j] << ' ' ;
    cout << endl ;
  }
}

void build()
{
  fr(i,1,n)
    fr(j,1,n){
    f[i][j]=2*(st[i][j]-st[i+1][j]);
    f[i][n+1]+=diff(st[i][j],st[i+1][j]);
  }
}

void solve()
{
  double mult;
  fr(i,1,n-1){
    fr(j,i+1,n){
      mult=f[j][i]/f[i][i];
      fr(p,i,n+1)
    f[j][p]-=mult*f[i][p];
    }
  }
  rf(i,1,n-1)
    fr(j,i+1,n){
    f[i][n+1]-=f[i][j]*f[j][n+1]/f[j][j];
    f[i][j]=0;
  }
}

int main()
{
#ifndef ONLINE_JUDGE
  freopen("1013.in","r",stdin);
  freopen("1013.out","w",stdout);
#endif
  cin >> n ;
  fr(i,1,n+1)
    fr(j,1,n)
    cin >> st[i][j] ;
  build();
  solve();
  fr(i,1,n-1)
    cout << fixed << setprecision(3) << f[i][n+1]/f[i][i] << '  ' ;
  cout << f[n][n+1]/f[n][n] ;
  return 0;
}