【高斯消元】BZOJ1013(JSOI2008)[球形空间产生器sphere]题解

题目概述

给出n+1个n维的点，求这几个点所在球的球心。

解题报告

一眼看去，好像没什么方法能够求出球心，那么只能用假设法了。设球心O的坐标为$(O_1,O_2,O_3,⋯,O_n)$，第i个点的坐标为$(x_{i,1},x_{i,2},x_{i,3},⋯,x_{i,n})$，根据欧几里得距离公式，就可以得出如下方程：
$∑(x_{1,i}^2-2x_{1,i} O_i+O_i^2)=∑(x_{2,i}^2-2x_{2,i}O_i+O_i^2)=⋯=∑(x_{n+1,i}^2-2x_{n+1,i} O_i+O_i^2)$
用1和2~n+1的等式得出n个方程：
$∑(2(x_{2,i}-x_{1,i})O_i )=∑(x_{2,i}^2-x_{1,i}^2)$
$∑(2(x_{3,i}-x_{1,i})O_i )=∑(x_{3,i}^2-x_{1,i}^2)$
$……$
$∑(2(x_{n+1,i}-x_{1,i})O_i )=∑(x_{n+1,i}^2-x_{1,i}^2)$
然后用高斯消元求解线性方程。

示例程序

#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
const int maxn=10;

int n,a[10],b[5];
double P[maxn+5][maxn+5],M[maxn+5][maxn+5],ans[maxn+5];

int fcmp(double x,double y) {if (fabs(x-y)<=1e-10) return -1;if (x<y) return 0; else return 1;}
int main()
{
    freopen("program.in","r",stdin);
    freopen("program.out","w",stdout);
    scanf("%d",&n);
    for (int i=1;i<=n+1;i++)
    for (int j=1;j<=n;j++)
        scanf("%lf",&P[i][j]);
    for (int i=2;i<=n+1;i++)
    for (int j=1;j<=n;j++)
        M[i-1][j]=(P[i][j]-P[1][j])*2,M[i-1][n+1]+=P[i][j]*P[i][j]-P[1][j]*P[1][j];
    //构造矩阵
    for (int t=1;t<=n;t++)
    {
        int where;
        for (int i=t;i<=n;i++) if (fcmp(M[i][t],0)!=-1) {where=i;break;} //找第t列不为0的行
        swap(M[where],M[t]); //与t行交换
        for (int i=t+1;i<=n;i++) //从t+1行开始，消元
        {
            double num=M[i][t]/M[t][t];
            for (int j=t;j<=n+1;j++) M[i][j]-=M[t][j]*num;
        }
    }
    for (int t=n;t>=1;t--)
    {
        ans[t]=M[t][n+1]; //取常数项
        for (int i=t+1;i<=n;i++) ans[t]-=M[t][i]*ans[i]; //代入
        ans[t]/=M[t][t];
    }
    for (int i=1;i<=n;i++) if (i==n) printf("%.3lf\n",ans[i]); else printf("%.3lf ",ans[i]);
    return 0;
}