[BZOJ 1013]球型空间产生器

不行不能颓了……线代其实刚起步……

注意每个点到球心的距离相同,距离之平方亦相同。然后假设所有半径平方为[tex]d[/tex],那么我们先可以列出两个方程(假设[tex]n=2[/tex],不过此时事实上可以列三个,但出于方便讨论的目的便只列两个):

[tex](x_1-a_1)^2+(x_2-a_2)^2=d[/tex]

[tex](x_1-b_1)^2+(x_2-b_2)^2=d[/tex]

两方程作差,得:

[tex](x_1-a_1)^2+(x_2-a_2)^2-(x_1-b_1)^2-(x_2-b_2)^2=0[/tex]

整理,得:

[tex]2(b_1-a_1)x_1+2(b_2-a_2)x_2=(b_1+a_1)(b_1-a_1)+(b_2+a_2)(b_2-a_2)[/tex]

对这种方法加以推广,便可列出[tex]n[/tex]个[tex]n[/tex]元一次方程。高斯消元求解即可。

代码:

/**************************************************************
    Problem: 1013
    User: danihao123
    Language: C++
    Result: Accepted
    Time:4 ms
    Memory:1296 kb
****************************************************************/
 
#include <cstdio>
#include <cmath>
#include <iostream>
#include <algorithm>
using namespace std;
const int maxn=20;
int n;
double M[maxn][maxn];
double A[maxn][maxn];
inline void Gause(){
    register int i,j,k,r;
    register double f;
    for(i=1;i<=n;i++){
        r=i;
        for(j=i+1;j<=n;j++)
            if(fabs(A[j][i])>fabs(A[r][i]))
                r=j;
        if(r!=i)
            for(j=1;j<=(n+1);j++)
                swap(A[r][j],A[i][j]);
        for(k=i+1;k<=n;k++){
            f=A[k][i]/A[i][i];
            for(j=i;j<=n+1;j++)
                A[k][j]-=f*A[i][j];
        }
    }
    for(i=n;i>=1;i--){
        for(j=i+1;j<=n;j++)
            A[i][n+1]-=A[j][n+1]*A[i][j];
        A[i][n+1]/=A[i][i];
    }
}
int main(){
    register int i,j;
    bool flag=false;
    cin>>n;
    for(i=1;i<=(n+1);i++)
        for(j=1;j<=n;j++)
            cin>>M[i][j];
    for(i=1;i<=n;i++)
        for(j=1;j<=(n+1);j++)
            A[i][j]=0;
    for(i=1;i<=n;i++){
        for(j=1;j<=n;j++){
            A[i][j]=2*(M[i+1][j]-M[i][j]);
            A[i][n+1]+=(M[i+1][j]-M[i][j])*(M[i+1][j]+M[i][j]);
        }
    }
    Gause();
    for(i=1;i<=n;i++){
        if(flag)
            putchar(' ');
        flag=true;
        printf("%.3f",A[i][n+1]);
    }
    putchar('\n');
    return 0;
}