[省选前题目整理][BZOJ 1013]球形空间产生器sphere(高斯消元)

题目链接

http://www.lydsy.com/JudgeOnline/problem.php?id=1013

思路

非常好的高斯消元练习题。。。
首先我们设在$n$维空间中，点a的坐标为$(a_1,a_2...a_n)$,点b的坐标为$(b_1,b_2...b_n)$，球心坐标$(x_1,x_2...x_n)$,球的半径为$r$
则
$(a_1-x_1)^2+(a_2-x_2)^2+...+(a_n-x_n)^2=r^2 ...①$
$(b_1-x_1)^2+(b_2-x_2)^2+...+(b_n-x_n)^2=r^2 ...②$
看上去这两个式子都是二次的，不能套用一次的高斯消元，但是当我们$①-②$后，神奇的一幕发生了：
$(a_1-b_1)x_1+(a_2-b_2)x_2+...+(a_n-b_n)x_n=\frac{(a_1^2-b_1^2)+(a_2^2-b_2^2)+...+(a_n^2-b_n^2)}{2}$
由于此题有$n$个未知数，我们需要$n$个方程。
以第一个点的坐标和第$2$~$n+1$个点相配对构造 n <script id="MathJax-Element-528" type="math/tex">n</script>个方程就ok了。常数项就是上面式子等号右边的部分
然后高斯消元搞定
蒟蒻我居然打错了高斯消元的交换行操作，WA了3发，我太弱了。。。

代码

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <algorithm>
#include <cmath>

#define MAXN 1100

using namespace std;

double pos[MAXN],a[MAXN][MAXN],ans[MAXN]; //pos[]保存第一个点的坐标，ans[]保存解

void Gauss(int n) //求n元一次方程
{
    for(int i=1;i<=n;i++) //上三角变换
    {
        int maxr=0;
        for(int j=i;j<=n;j++)
            if(fabs(a[j][i])>fabs(a[maxr][i]))
                maxr=j;
        for(int j=1;j<=n+1;j++) //!!!!!交换行，使得第i行的第i个元素是所有行里最大的
            swap(a[i][j],a[maxr][j]);
        for(int j=i+1;j<=n;j++) //将第j行的第i个元素消去
        {
            double tmp=-a[j][i]/a[i][i];
            for(int k=i;k<=n+1;k++) //!!!!!!!
                a[j][k]+=a[i][k]*tmp;
        }
    }
    for(int i=n;i>=1;i--) //回代
    {
        for(int j=n;j>i;j--) //第i行的i+1~n列都要消去
            a[i][n+1]-=a[i][j]*ans[j];
        ans[i]=a[i][n+1]/a[i][i]; //!!!!!!
    }
}

int main()
{
    int n;
    scanf("%d",&n);
    for(int i=1;i<=n;i++) scanf("%lf",&pos[i]);
    for(int i=1;i<=n;i++)
    {
        double tmp[MAXN];
        for(int j=1;j<=n;j++)
            scanf("%lf",&tmp[j]);
        for(int j=1;j<=n;j++)
        {
            a[i][j]=pos[j]-tmp[j];
            a[i][n+1]+=(pos[j]*pos[j]-tmp[j]*tmp[j])/2;
        }
    }
    Gauss(n);
    for(int i=1;i<=n;i++)
        printf("%.3lf%c",ans[i],i==n?'\n':' ');
    return 0;
}