BZOJ 1013 高斯消元

对于二维空间中的圆，有 ((x-a)^2 + (y-b)^2 = r^2)
对于三维空间中的球体，有 ((x-a)^2 + (y-b)^2 + (z - c)^2 = r^2)
设球心的坐标为((x_1, x_2, …, x_n))
对于(n)维空间中的球体，对于每个(j)，有(Sigma_{i}(x_i - a_{ji})^2 = r^2)。设它为第(j)个方程。
可得到(n)个未知数，(n+1)个二次方程
将第(j)个方程和第(j-1)个方程相减，得(n)个一次方程。
套上高斯消元版子就行咯。

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
  
 #define rep(i,a,b) for (int i = a; i <= b; i++)
 #define dep(i,a,b) for (int i = a; i >= b; i--) 
  
 const int N = 15;
  
 int n;
 double s[N], a[N][N];
  
 void init() {
    double t;
    rep(j,1,n) scanf("%lf", &s[j]), a[j][n+1] = 0;
    rep(i,1,n) 
        rep(j,1,n) { 
            scanf("%lf", &t);
            a[i][j] = 2*(t - s[j]);
            a[i][n+1] += -(s[j]*s[j] - t*t);
        }
 }
  
 void Gauss() {
    int now = 1;
    double t;
    rep(i, 1, n) {
        now = i;
        rep(j, i, n) if (a[j][i] > a[now][i]) now = j;
        rep(j, 1, n+1) swap(a[now][j], a[i][j]);
        t = a[i][i];
        rep(j, i, n+1) a[i][j] /= t;
        rep(j, i+1, n) {
            t = a[j][i];
            rep(k, i, n+1) a[j][k] -= t*a[i][k];
        }
    }
    dep(i, n, 1) rep(j, i+1, n) a[i][n+1] -= a[j][n+1]*a[i][j];
 }
  
int main()
{
    scanf("%d", &n);
    init();
    Gauss();
  
    rep(i,1,n-1) printf("%.3f ", a[i][n+1]);
    printf("%.3f
", a[n][n+1]);
  
    return 0;
}