HYSBZ/BZOJ 1013 [JSOI2008] 球形空间产生器sphere - 高斯约当消元

来源:互联网 发布:国动网络通信集团图片 编辑:程序博客网 时间:2024/06/04 18:55

题目描述

分析:

根据题目给出的dist的定义,并设球心的坐标为(x1,x2,...xn),列出n+1个方程。
step 1. n+1个方程的两边同时平方。
step 2. n+1个方程拆括号,化简。
step 3. 前n个方程-第n+1个方程,得到新的n个方程(dist2)
step 4. n个关于n个未知数的方程,解吧。

#include<cstdio>#include<cmath>#include<algorithm>using namespace std;#define MAXN 10const double eps=1e-12;int n;double a[MAXN+10][MAXN+10],b[MAXN+10][MAXN+10];void read(){    scanf("%d",&n);    for(int i=1;i<=n+1;i++){        for(int j=1;j<=n;j++)            scanf("%lf",&b[i][j]);    }    for(int i=1;i<=n;i++){        a[i][n+1]=0.0;        for(int j=1;j<=n;j++){            a[i][j]=2.0*b[n+1][j]-2.0*b[i][j];            a[i][n+1]+=b[n+1][j]*b[n+1][j]-b[i][j]*b[i][j];        }    }}void GJ_elimination(int equ,int var,int &row,int &col){    int mx;    for(row=col=1;row<=equ&&col<=var;row++,col++){        mx=row;        for(int i=row+1;i<=equ;i++)            if(fabs(a[i][col])>fabs(a[mx][col]))                mx=i;        if(mx!=row)            swap(a[mx],a[row]);        if(fabs(a[row][col])<eps){            row--;            continue;        }        for(int i=1;i<=equ;i++){            if(i==row||fabs(a[i][col])<eps)                continue;            for(int j=var+1;j>=col;j--)                a[i][j]-=a[i][col]/a[row][col]*a[row][j];        }    }}void Gauss_Jordan(int equ,int var){    int row,col;    GJ_elimination(equ,var,row,col);    for(int i=1;i<equ;i++)        printf("%.3lf ",a[i][var+1]/a[i][i]);    printf("%.3lf\n",a[equ][var+1]/a[equ][equ]);}int main(){    read();    Gauss_Jordan(n,n);}
0 0
原创粉丝点击