[BZOJ 1013][JSOI 2008] 球形空間產生器sphere 題解（高斯消元）

[BZOJ 1013][JSOI 2008] 球形空間產生器sphere

Description

有一個球形空間產生器可以在n維空間中產生一個堅硬的球體。如今，你被困在了這個n維球體中，你只知道球
面上n+1個點的座標，你須要以最快的速度肯定這個n維球體的球心座標，以便於摧毀這個球形空間產生器。

Input
第一行是一個整數n(1<=N=10)。接下來的n+1行，每行有n個實數，表示球面上一點的n維座標。每個實數精確到小數點
後6位，且其絕對值都不超過20000。

Output
有且只有一行，依次給出球心的n維座標（n個實數），兩個實數之間用一個空格隔開。每一個實數精確到小數點
後3位。數據保證有解。你的答案必須和標準輸出如出一轍纔可以得分。

Solution

1.考慮構造對於不一樣的店一樣結構的方程。

由於是一個球，因此每一個點到球心的距離都相等，

咱們設這個半徑爲R，球心座標爲O(x1,x2,....,xn);

那麼對於每個點P(ai1,ai2,...,ain)：咱們易得

sqrt ( ( ai1 - x1 ) ^ 2 + ( ai2 - x2 ) ^ 2 + ... + ( ain - xn ) ^2 ) = R;

2.考慮構造方程組

將上式兩側平方再展開，得

-2 ( ai1 * x1 + ai2 * x2 + ... + ain * xn )
+( ai1 ^ 2 + ai2 ^ 2 + ... + ain ^ 2 )
+( x1 ^ 2 + x2 ^ 2 + ... + xn ^ 2 )
= R ^ 2;

這時咱們看數據，給出n+1個點，n個點就能夠構造該方程組，那多給的一個點是用來幹什麼的呢（設選第一個點爲這個點）？

沒錯！消掉重複出現的部分( x1 ^ 2 + x2 ^ 2 + ... + xn ^ 2 ) 和R ^ 2，

即令其餘全部的點的方程都減掉多的一個點的方程，整理獲得其餘方程格式爲：

2 ( ai1 * x1 + ai2 * x2 + ... + ain * xn )
= ( ai1 ^ 2 + ai2 ^ 2 + ... + ain ^ 2 )
- ( a11 ^ 2 + a12 ^ 2 + ... + a1n ^ 2 )
- 2 ( a11 * x1 + a12 * x2 + ... + a1n * xn ) ;

右側是常數，左側展開就是一個愉快的高斯消元方程組，解方程組便可。

這裏使用的高斯消元法是一種比較毒瘤的方法，詳解參考個人隨筆：http://www.cnblogs.com/COLIN-LIGHTNING/p/8981923.html

Code

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
using namespace std;

const int max_n=15;
double a[max_n][max_n+1],v[max_n],del[max_n],x;
int n,w[max_n]; 

inline double sqr(double x){return x*x;} 

void init(){
    for(int i=1;i<=n;++i)scanf("%lf",&del[i]);
    for(int i=1;i<=n;++i){
        for(int j=1;j<=n;++j){
            scanf("%lf",&x);
            a[i][n+1]+=sqr(x)-sqr(del[j]);
            a[i][j]=2*(x-del[j]);
        }
    }
}

void gauss(){
    double eps=1e-6;
    for(int i=1;i<=n;++i){//enumerate the equation;
        int p=0;          //Record the position of the largest number; 
        double mx=0;      //Recording the largest number;
        for(int j=1;j<=n;++j)
            if(fabs(a[i][j])-eps>mx){
                mx=fabs(a[i][j]);p=j;//fabs() returns the absolute value of float; 
            }
        w[i]=p;
        for(int j=1;j<=n;++j)
            if(i!=j){       //other equations
                double t=a[j][p]/a[i][p];
                for(int k=1;k<=n+1;++k)//n+1 is important
                    a[j][k]-=a[i][k]*t;
            }
    }
    for(int i=1;i<=n;++i) v[w[i]]=a[i][n+1]/a[i][w[i]];
    for(int i=1;i<=n;++i) printf("%.3lf ",v[i]);
}

int main(){
    scanf("%d",&n);
    init();
    gauss();
    return 0; 
}