洛谷P4035 [JSOI2008]球形空间产生器(高斯消元)

洛谷题目传送门

球啊球 @xzz_233 qaq

高斯消元模板题,关键在于将已知条件转化为方程组。

可以发现题目要求的未知量有(n)个,题目却给了我们(n+1)个点的坐标,这其中必有玄机。

由高中数学知识可以知道,三点定圆(二维),四点定球(三维)······以此类推,应该是(n+1)个点才能确定一个(n)维空间下的球。

那么隐藏的另一个关键未知量在哪里呢?

想想圆的标准方程((x-x_0)^2+(y-y_0)^2=r^2),除了圆心坐标,半径不也对这个圆起到决定性作用么?

接下来,额外设一个未知量——球的半径(r),开始试着对条件式进行变换。

对于(n+1)个点,它们与球心的距离是定值(r),那么我们可以得到形式如下的(n+1)个方程(a为球面点坐标,x为球心坐标)

[(a_1-x_1)^2+(a_2-x_2)^2+···+(a_n-x_n)^2=r^2 ]

显然我们要把已知量和未知量分开,于是展开,移项

[a_1^2-2a_1x_1+x_1^2+a_2^2-2a_2x_2+x_2^2+···+a_n^2-2a_nx_n+x_n^2=r^2 ]

[2a_1x_1+2a_2x_2+···+2a_nx_n+r^2-x_1^2-x_2^2-···-x_n^2=a_1^2+a_2^2+···+a_n^2 ]

发现(r^2-x_1^2-x_2^2-···-x_n^2)(a)无关,所以考虑换元,设(t=r^2-x_1^2-x_2^2-···-x_n^2)(实际上我们并不用求(r)

终于,我们可以看到一个关于(x_1,x_2,···,x_n,t)((n+1))元方程组了,上高斯消元

具体实现看代码

#include<cmath>
#include<cstdio>
#define R register
#define init for(i=1;i<=n;++i)ne[i-1]=pr[i+1]=i
const int N=19;
int p[N],pr[N],ne[N];
double a[N][N];
int main(){
    R int n,i,j,k,x;
    R double mx,d;
    scanf("%d",&n);++n;
    init;//链表初始化,为了实现交换行
    for(i=1;i<=n;++i){
        for(j=1;j<n;++j){
            scanf("%lf",&d);//处理系数
            a[i][j]=d*2;a[i][n+1]+=d*d;
        }
        a[i][n]=1;//t的系数为1
    }
    for(k=1;k<=n;++k){
        mx=0;//蒟蒻没有交换主元,而是交换行
                //这样做防掉精度的效果可能不如交换主元
        for(i=ne[0];i;i=ne[i])
            if(mx<fabs(a[i][k]))
                mx=fabs(a[i][k]),x=i;
        d=a[p[k]=x][k];//选择当前a最大的一行
        pr[ne[pr[x]]=ne[x]]=pr[x];
        for(j=1;j<=n+1;++j)
            a[x][j]/=d;
        for(i=ne[0];i;i=ne[i])
            for(d=a[i][j=k];j<=n+1;++j)
                a[i][j]-=d*a[x][j];
    }//高斯消元
    init;
    for(k=n;k;--k){
        d=a[x=p[k]][n+1];
        pr[ne[pr[x]]=ne[x]]=pr[x];
        for(i=ne[0];i;i=ne[i])
            a[i][n+1]-=d*a[i][k];
    }//回代
    for(k=1;k<n;++k)
        printf("%.3f ",a[p[k]][n+1]);
    puts("");
    return 0;
}
原文地址:https://www.cnblogs.com/flashhu/p/9193318.html