BZOJ 1013 球形空间产生器sphere (double高斯消元)
1013: [JSOI2008]球形空间产生器sphere
Time Limit: 1 Sec Memory Limit: 162 MBSubmit: 2696 Solved: 1397
[Submit][Status]
Description
有一个球形空间产生器能够在n维空间中产生一个坚硬的球体。现在，你被困在了这个n维球体中，你只知道球面上n+1个点的坐标，你需要以最快的速度确定这个n维球体的球心坐标，以便于摧毁这个球形空间产生器。
Input
第一行是一个整数，n。接下来的n+1行，每行有n个实数，表示球面上一点的n维坐标。每一个实数精确到小数点后6位，且其绝对值都不超过20000。
Output
有且只有一行，依次给出球心的n维坐标（n个实数），两个实数之间用一个空格隔开。每个实数精确到小数点后3位。数据保证有解。你的答案必须和标准输出一模一样才能够得分。
Sample Input
0.0 0.0
-1.0 1.0
1.0 0.0
Sample Output
HINT
数据规模：
对于40%的数据，1<=n<=3
对于100%的数据，1<=n<=10
提示：给出两个定义：
1、 球心：到球面上任意一点距离都相等的点。
2、 距离：设两个n为空间上的点A, B的坐标为(a1, a2, …, an), (b1, b2, …, bn)，则AB的距离定义为：dist = sqrt( (a1-b1)^2 + (a2-b2)^2 + … + (an-bn)^2 )
传送门：http://www.lydsy.com/JudgeOnline/problem.php?id=1013
思路：先从二维进行考虑, 设圆心（x，y），给定的点（a，b）,（a，b）到圆心的距离为 (a-x)^2+(b-y)^2 = a^2+2ax+x^2+b^2+2by+y^2, 于是我们可以用一个点将其它两个点变为俩个方程,例如还有一个点（a1，b1）,则2(a1-a)x+2(b1-b)y=a1^2-a^2+b1^2-b^2,然后解方程。
我本来想用bin神的double 高斯消元模版的...但我好象不会用...结果一直都为0...希望有人能帮我看下我哪里错了...
代码1：(用bin神板子..有错！)
#include<algorithm>
#include<iostream>
#include<string.h>
#include<math.h>
using namespace std;
#define eps 1e-9
const int MAXN=220;
double a[MAXN][MAXN],x[MAXN];//方程的左边的矩阵和等式右边的值，求解之后x存的就是结果
void Debug(int equ,int var)
{
int i, j;
for (i = 0; i < equ; i++)
{
for (j = 0; j < var + 1; j++)
{
cout << a[i][j] << " ";
}
cout << endl;
}
cout << endl;
}
/*
*返回0表示无解，1表示有解
*/
int Gauss(int equ,int var)//方程数和未知数个数
{
int i,j,k,col,max_r;
for(k=0,col=0; k<equ&&col<var; k++,col++)
{
max_r=k;
for(i=k+1; i<equ; i++)
if(fabs(a[i][col])>fabs(a[max_r][col]))
max_r=i;
if(fabs(a[max_r][col])<eps)return 0;
if(k!=max_r)
{
for(j=col; j<var; j++)
swap(a[k][j],a[max_r][j]);
swap(x[k],x[max_r]);
}
x[k]/=a[k][col];
for(j=col+1; j<var; j++)a[k][j]/=a[k][col];
a[k][col]=1;
for(i=0; i<equ; i++)
if(i!=k)
{
x[i]-=x[k]*a[i][k];
for(j=col+1; j<var; j++)a[i][j]-=a[k][j]*a[i][col];
a[i][col]=0;
}
}
return 1;
}
double b[MAXN][MAXN];
int main()
{
freopen("in.txt", "r", stdin);
int equ,var,n;
scanf("%d",&n);
equ = var = n;
for(int i = 0; i <= n; i++)
for(int j = 0; j < n; j++)
scanf("%lf",&b[i][j]);
memset(a, 0, sizeof(a));
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
a[i][j] = 2*(b[i+1][j]-b[i][j]);
a[i][n] += b[i+1][j]*b[i+1][j] - b[i][j]*b[i][j];
}
}
//Debug(equ,var);
Gauss(equ,var);
for (int i = 0; i < var-1; i++)
printf("%.3lf ",x[i]);
printf("%.3lf\n",x[var-1]);
return 0;
}
代码2：(用别人的板子)
#include<algorithm>
#include<iostream>
#include<string.h>
#include<math.h>
using namespace std;
#define eps 1e-8
const int MAXN=220;
double f[MAXN][MAXN], x[MAXN];
void Gauss(int n)
{
for (int i = 1; i <= n - 1; ++i)
{
int maxi = i;
for (int j = i + 1; j <= n; ++j)
if (fabs(f[j][i]) > fabs(f[maxi][i])) maxi = j;
if (fabs(f[maxi][i]) > eps)
{
for (int j = 1; j <= n + 1; ++j) swap(f[maxi][j], f[i][j]);
for (int j = i + 1; j <= n; ++j)
if (fabs(f[j][i]) > eps)
{
double gs = -f[j][i] / f[i][i];
for (int k = i; k <= n + 1; ++k)
f[j][k] += f[i][k] * gs;
}
}
}
for (int i = n; i >= 1; --i)
{
double res = f[i][n + 1];
for (int j = i + 1; j <= n; ++j)
res -= x[j] * f[i][j];
x[i] = res / f[i][i];
}
}
double b[MAXN][MAXN];
int main()
{
//freopen("in.txt", "r", stdin);
int n;
scanf("%d",&n);
for(int i = 1; i <= n+1; i++)
for(int j = 1; j <= n; j++)
scanf("%lf",&b[i][j]);
memset(f, 0, sizeof(f));
for (int i = 1; i <= n; i++)
{
for (int j = 1; j <= n; j++)
{
f[i][j] = 2*(b[i+1][j]-b[i][j]);
f[i][n+1] += b[i+1][j]*b[i+1][j] - b[i][j]*b[i][j];
}
}
Gauss(n);
for (int i = 1; i <= n - 1; ++i)
printf("%.3lf ", x[i]);
printf("%.3lf", x[n]);
return 0;
}
