【BZOJ】1013: [JSOI2008]球形空间产生器sphere

【BZOJ】1013: [JSOI2008]球形空间产生器sphere

题意:给n+1个n维的点的坐标，要你求出一个到这n+1个点距离相等的点的坐标；

思路:高斯消元即第i个点和第i+1个点处理出一个式子，这样n+1个点正好有n个系数的n元变量，即可求解。

式子:Σ( (a[i][j] - x[j])^2 )  = Σ( a[i+1][j] - x[j])^2 )

=>   Σ( x[j]*[2*(a[i+1][j]-a[i][j])] ) = Σ（a[i+1][j]*a[i+1][j] - a[i][j]*a[i][j]);直接预处理即可；

注意:在Gauss处理出上三角阵的过程中，每次要选出主对角线绝对值最大的行作为参考行，貌似是精度问题。还有就是归零的过程中，要变成参考行再消，为了不出现除0的情况。

#include<iostream>

#include<cstdio>

#include<cstring>

#include<string.h>

#include<algorithm>

#include<map>

#include<queue>

#include<vector>

#include<cmath>

#include<stdlib.h>

#include<time.h>

using namespace std;

typedef long long ll;

#define rep0(i,l,r) for(int i = (l);i < (r);i++)

#define rep1(i,l,r) for(int i = (l);i <= (r);i++)

#define MS0(a) memset(a,0,sizeof(a))

#define MS1(a) memset(a,-1,sizeof(a))

double a[][],A[][];

int n;

void Gauss()

{

    int i,j,k;

    rep1(i,,n){

        int mx = i;

        rep1(j,i+,n) if(fabs(A[mx][i]) < fabs(A[j][i])) mx = j;

        rep1(j,i,n+) swap(A[mx][j],A[i][j]);

        rep1(j,i+,n)if(A[i][i] != ){

            double y = A[j][i]/A[i][i];

            rep1(k,i,n+) A[j][k] -= y*A[i][k];

        }

    }

    for(int i = n;i >= ;i--){

        rep1(j,i+,n)   A[i][n+] -= A[i][j] * A[j][n+];

        A[i][n+] /= A[i][i];   //化为系数为1；保证有解，则A[i][i] != 0;

    }

}

int main()

{

    int i,j;

    scanf("%d",&n);

    rep1(i,,n+)

        rep1(j,,n)

            scanf("%lf",&a[i][j]);

    rep1(i,,n)

        rep1(j,,n){

            A[i][j] = *(a[i+][j] - a[i][j]);

            A[i][n+] += a[i+][j]*a[i+][j] - a[i][j]*a[i][j];

        }

    Gauss();

    printf("%.3f",A[][n+]);

    rep1(i,,n) printf(" %.3f",A[i][n+]);

}