bzoj1013 [JSOI2008]球形空间产生器sphere (高斯消元)

1013  [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

的坐标为(a1, a2, …, an), (b1, b2, …, bn)，则AB的距离定义为：dist = sqrt( (a1-b1)^2 + (a2-b2)^2 + 

… + (an-bn)^2 )

---

传送门

 

事情是这样的，本来是在写另外一题，肝了半天没写出来，查了下题号发现标签是“概率DP+高斯消元”。由于姿势水平不够，并不知道这是个什么骚东西，找道题来学一下。

 

（UPD：最后还是没写出来。。。

 

高斯消元是个用于解n元线性方程组的玩意。

对于线性方程组$\begin{cases} a_1 x_1 + a_2 x_2 +a_3 x_3 = d_1\\b_1 x_1 + b_2 x_2 +b_3 x_3 = d_2\\c_1 x_1 + c_2 x_2 +c_3 x_3 = d_3 \end{cases}$

我们可以写出一个等价的系数矩阵$\begin{bmatrix}a_1&a_2&a_3\\b_1&b_2&b_3\\c_1&c_2&c_3\end{bmatrix} \times \begin{bmatrix}x_1\\x_2\\x_3\end{bmatrix} = \begin{bmatrix}d_1\\d_2\\d_3\end{bmatrix}$

而我们希望得到的是这样一个矩阵：$\begin{bmatrix}a_{11}&a_{12}&a_{13}\\0&a_{22}&a_{23}\\0&0&a_{33}\end{bmatrix} \times \begin{bmatrix}x_1\\x_2\\x_3\end{bmatrix} = \begin{bmatrix}b_1\\b_2\\b_3\end{bmatrix}$

系数矩阵的最后一行是$x_3$的解，只要逐层向上代入消元就能解出所有未知数。

具体的算法流程是这样的：

1 枚举矩阵中的元素$a_{i,i}$，

　　1.1 从其下方中找到绝对值最大的元素$a_{j,i}$；

　　1.2 将第$i$行与第$j$行$[i,n+1]$的部分交换；

　　1.3 将第$i$列第$[i+1,n]$行调整为$0$；

2 从第$n$行至第$1$行进行消元；

 

对于这道题，距离公式带平方，但是题目给了$n+1$条式子，用后$n$条式子减去第一条式子即可消去二次项。

#include<cstring>
#include<cstdio>
#include<cmath>
#include<algorithm>
#define foru(i,x,y) for(int i=x;i<=y;i++)
#define ford(i,x,y) for(int i=x;i>=y;i--)
using namespace std;
double f[100][100],x,r[100];
int n;

void gauss(){
    foru(i,1,n){
        int t=i;
        foru(j,i+1,n)if(fabs(f[j][i])>fabs(f[t][i]))t=j;
        if(t!=i)foru(j,i,n+1)swap(f[i][j],f[t][j]);//交换 提高精度 
        foru(j,i+1,n){
            double x=f[j][i]/f[i][i];
            foru(k,i,n+1)f[j][k]-=x*f[i][k];
            //f[i][i]以下的元素用f[i][i]按比例减为0,同行的其它元素按比例减去第一行的对应值 
            //即矩阵的初等行变换 
        }
    }
    ford(i,n,1){
        double t=(f[i][n+1]/=f[i][i]);
        foru(j,1,i-1)f[j][n+1]-=(t*f[j][i]);//消元 
    }
}

int main(){
    scanf("%d",&n);
    foru(i,1,n)scanf("%lf",&r[i]);
    foru(i,1,n){
        foru(j,1,n){
            scanf("%lf",&x);
            f[i][j]=2*(x-r[j]);
            f[i][n+1]+=x*x-r[j]*r[j];//消去二次项 
        }
    }
    gauss();
    printf("%.3lf",f[1][n+1]);
    foru(i,2,n)printf(" %.3lf",f[i][n+1]);
    return 0; 
}