P4035

· · 个人记录

[JSOI2008]球形空间产生器

高斯消元。

简单理解就是将线性方程组的系数化成矩阵，再与右侧常数项共同构成一个增广矩阵。

然后就是机械消元，通过一个方程乘上某个非零系数把其余方程的此项全部化 0，然后看下一项再重复执行次操作。

这题不存在主元和自由元，所以解出来就完事了。

这题所谓圆心的含义其实就是所有点到该点的距离相同。

所以就会有方程（组）：

\sum_{j=1}^n(a_i,_j-x_j)^2=C

其中 C 为定值，是个常数。

然后这里我们尝试把这 n+1 个非线性方程转化为线性方程组，显然作差是一个不错的选择。

那么相邻的两个方程作差，我们就可以得到含有 n 个方程的线性方程组：

\sum_{j=1}^n(a^2_i,_j+x_j^2-2x_j\cdot a_i,_j)-\sum_{j=1}^n(a^2_{i+1},_j+x_j^2-2x_j\cdot a_{i+1},_j)=0 \sum_{j=1}^n(a^2_i,_j-a^2_{i+1},_j-2x_j(a_i,_j-a_{i+1},_j))=0 \sum_{j=1}^n2(a_i,_j-a_{i+1},_j)x_j=\sum_{j=1}^n(a^2_i,_j-a^2_{i+1},_j)

这里 i=1,2,\cdots,n。

然后我们就转化成下面的增广矩阵：

\begin{bmatrix}2(a_1,1-a_2,_1)&2(a_1,_2-a_2,_2)&\cdots&2(a_1,_n-a_2,n)&\sum_{j=1}^n(a^2_1,_j-a^2_2,_j)\\2(a_2,1-a_3,_1)&2(a_2,_2-a_3,_2)&\cdots&2(a_2,_n-a_3,n)&\sum_{j=1}^n(a^2_2,_j-a^2_3,_j)\\\vdots&\vdots&\ddots&\vdots&\vdots\\2(a_n,1-a_{n+1},_1)&2(a_n,_2-a_{n+1},_2)&\cdots&2(a_n,_n-a_{n+1},n)&\sum_{j=1}^n(a^2_n,_j-a^2_{n+1},_j)\end{bmatrix}

然后消元就完事了。

时间复杂度 O(n^3)。

代码：

#include<iostream>
#include<cstdio>
#include<cmath>
#define ll long long
using namespace std;

const ll N=20;

const double acc=1e-8;

ll n;

double a[N+5][N+5],b[N+5],c[N+5][N+5];

int main() {

    scanf("%lld",&n);

    for(ll i=1;i<=n+1;i++) {
        for(ll j=1;j<=n;j++) {
            scanf("%lf",&a[i][j]);
        }
    }

    for(ll i=1;i<=n;i++) {
        for(ll j=1;j<=n;j++) {
            c[i][j]=2*(a[i][j]-a[i+1][j]);
            b[i]+=a[i][j]*a[i][j]-a[i+1][j]*a[i+1][j];
        }
    }

    for(ll i=1;i<=n;i++) {
        for(ll j=i;j<=n;j++) {
            if(fabs(c[j][i])>acc) {
                for(ll k=1;k<=n;k++) swap(c[i][k],c[j][k]);
                swap(b[i],b[j]);
            }
        }
        for(ll j=1;j<=n;j++) {
            if(i==j) continue;
            double rate=c[j][i]/c[i][i];
            for(ll k=i;k<=n;k++) c[j][k]-=c[i][k]*rate;
            b[j]-=b[i]*rate;
        }
    }

    for(ll i=1;i<=n;i++) printf("%.3f ",b[i]/c[i][i]);

    return 0;
}