BZOJ - 1013 高斯消元

Posted 2020-10-22 The Azure Arbitrator

n维空间中给出n+1个球面上的点求圆心坐标(x0,x1,...xn-1)
任选其中一个点坐标如第一个点(a0,b0...z0)
(x0-a0)^2+(x1-b0)^2+...=r^2
对于剩下的n个点都与上面的式子作差,把高次方的消去得到线性方程组
2(a1-a0)x0+2(b1-b0)x1+2(c1-c0)x2+...=dis1-dis0
...
共n行n列,
然后构造Ax=b丢进去高斯消元就行了
//其实我只是来试试板子(来自挑战程序设计竞赛)
//用vector来构造矩阵有点不适应

#include<bits/stdc++.h>
#define rep(i,j,k) for(int i=j;i<=k;i++)
using namespace std;
const int maxn = 1e5+11;
const double eps = 1e-8;
typedef vector<double> vec;
typedef vector<vec> mat;
vec gauss(const mat& A,const vec& b){
    int n=A.size();
    mat B(n,vec(n+1));//n element with vec(n+1)
    rep(i,0,n-1)rep(j,0,n-1)B[i][j]=A[i][j];
    rep(i,0,n-1) B[i][n]=b[i];//增广 
    rep(i,0,n-1){//未知数系数绝对值最大的移至第i行 
        int now=i;
        rep(j,i,n-1){
            if(abs(B[j][i])>abs(B[now][i])){
                now=j;
            }
        }
        swap(B[i],B[now]);
        if(abs(B[i][i])<eps)return vec();
        rep(j,i+1,n) B[i][j]/=B[i][i];
        rep(j,0,n-1) if(i!=j){
            rep(k,i+1,n){
                B[j][k]-=B[j][i]*B[i][k];
            }
        } 
    }
    vec x(n);
    rep(i,0,n-1) x[i]=B[i][n];
    return x; 
}
int main(){
    int n;
    double a0[22];
    while(scanf("%d",&n)!=EOF){
        mat A(n,vec(n));//n行n列 
        vec b(n);//n行1列 
        rep(i,0,n-1) scanf("%lf",&a0[i]);
        double dis0=0;
        rep(i,0,n-1) dis0+=a0[i]*a0[i];
        rep(i,0,n-1){
            double dis2=0;
            rep(j,0,n-1){
                scanf("%lf",&A[i][j]);
                dis2+=A[i][j]*A[i][j];
                A[i][j]=2*(A[i][j]-a0[j]);
            }
            b[i]=dis2-dis0;
        }
        vec x=gauss(A,b);
        for(int i=0;i<x.size();i++){
            if(i==x.size()-1) printf("%.3lf\n",x[i]);
            else printf("%.3lf ",x[i]);
        }
    }
    return 0;
}