BZOJ 1013: [JSOI2008]球形空间产生器sphere 高斯消元

原题链接：http://www.lydsy.com/JudgeOnline/problem.php?id=1013

题解：

由所有点到圆心的距离都是相等的，我们可以列出n个方程式，其中第$i$个是这样的：

$$\sum_{j=1}^{n} 2*(a_{vj}-a_{uj})*P_j=\sum_{j=1}^{n}-a_{uj}^2+a_{vj}^2$$

其中$v=i+1,u=i$

我们得到的解是$P$矩阵

代码：

/**************************************************************
Problem: 1013
User: HarryGuo2012
Language: C++
Result: Accepted
Time:0 ms
Memory:1276 kb
****************************************************************/

#include<iostream>
#include<vector>
#include<cmath>
#include<cstring>
#include<iomanip>
#define MAX_N 15
using namespace std;

const double EPS=1e-8;

double a[MAX_N][MAX_N];

int n;

const int maxn=MAX_N;
const double eps=1e-12;
int equ,var;//equ个方程,var个变量
double x[maxn];//解集
bool free_x[maxn];

int sgn(double x)
{
return (x>eps)-(x<-eps);
}
// 高斯消元法解方程组(Gauss-Jordan elimination).(0表示无解，1表示唯一解，大于1表示无穷解，并返回自由变元的个数)
int gauss() {
equ = n, var = n;//多少个方程，多少个变量
int i, j, k;
int max_r; // 当前这列绝对值最大的行.
int col; // 当前处理的列.
double temp;
int free_x_num;
int free_index;
// 转换为阶梯阵.
col = 0; // 当前处理的列.
memset(free_x, true, sizeof(free_x));
for (k = 0; k < equ && col < var; k++, col++) {
max_r = k;
for (i = k + 1; i < equ; i++) {
if (sgn(fabs(a[i][col]) - fabs(a[max_r][col])) > 0)
max_r = i;
}
if (max_r != k) { // 与第k行交换.
for (j = k; j < var + 1; j++)
swap(a[k][j], a[max_r][j]);
}
if (sgn(a[k][col]) == 0) { // 说明该col列第k行以下全是0了，则处理当前行的下一列.
k--;
continue;
}
for (i = k + 1; i < equ; i++) { // 枚举要删去的行.
if (sgn(a[i][col]) != 0) {
temp = a[i][col] / a[k][col];
for (j = col; j < var + 1; j++) {
a[i][j] = a[i][j] - a[k][j] * temp;
}
}
}
}

for (i = k; i < equ; i++) {
if (sgn(a[i][col]) != 0)
return 0;
}
if (k < var) {
for (i = k - 1; i >= 0; i--) {
free_x_num = 0;
for (j = 0; j < var; j++) {
if (sgn(a[i][j]) != 0 && free_x[j])
free_x_num++, free_index = j;
}
if (free_x_num > 1) continue;
temp = a[i][var];
for (j = 0; j < var; j++) {
if (sgn(a[i][j]) != 0 && j != free_index)
temp -= a[i][j] * x[j];
}
x[free_index] = temp / a[i][free_index];
free_x[free_index] = 0;
}
return var - k;
}

for (i = var - 1; i >= 0; i--) {
temp = a[i][var];
for (j = i + 1; j < var; j++) {
if (sgn(a[i][j]) != 0)
temp -= a[i][j] * x[j];
}
x[i] = temp / a[i][i];
}
return 1;
}

double point[MAX_N][MAX_N];

int main() {
cin >> n;
for (int i = 0; i < n + 1; i++)
for (int j = 0; j < n; j++)
cin >> point[i][j];
for (int i = 0; i < n; i++) {
double tmp = 0;
for (int j = 0; j < n; j++) {
a[i][j] = 2 * (point[i + 1][j] - point[i][j]);
tmp += point[i + 1][j] * point[i + 1][j] - point[i][j] * point[i][j];
}
a[i]
= tmp;
}
gauss();
for (int i = 0; i < n; i++) {
cout << setprecision(3) << fixed << x[i];
if (i != n - 1)cout << " ";
}
cout << endl;
return 0;
}