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][n] = tmp;
    }
    gauss();
    for (int i = 0; i < n; i++) {
        cout << setprecision(3) << fixed << x[i];
        if (i != n - 1)cout << " ";
    }
    cout << endl;
    return 0;
}