球形空间产生器

AcWing-207-球形空间产生器

分析

见《进阶指南》第161页。

实现

#include <iostream>
#include <cmath>
using namespace std;

const int N = 20;
const double epx = 1e-8;
int n;
double a[N][N], b[N][N], d[N]; // 增广矩阵：[b, d^T]

void Gauss () { // 数据保证有唯一解
    for (int c = 1, r = 1; c <= n; ++ c, ++ r) {
        int idx = r;
        for (int i = r; i <= n; ++ i)
            if (fabs(b[i][c]) > fabs(b[idx][c]))
                idx = i;
        swap(b[r], b[idx]);
        swap(d[r], d[idx]);

        for (int i = 1; i <= n; ++ i) {
            if (i == r) continue;
            double rate = b[i][c] / b[r][c];
            for (int j = c; j <= n; ++ j)
                b[i][j] -= rate * b[r][j];
            d[i] -= rate * d[r];
        }
    }
    for (int i = 1; i <= n; ++ i) d[i] /= b[i][i];
}

int main () {
    scanf("%d", &n);
    for (int i = 1; i <= n + 1; ++ i)
        for (int j = 1; j <= n; ++ j)
            scanf("%lf", &a[i][j]);
    
    for (int i = 1; i <= n; ++ i) {
        for (int j = 1; j <= n; ++ j) {
            b[i][j] = 2 * (a[i][j] - a[i + 1][j]);
            d[i] += a[i][j] * a[i][j] - a[i + 1][j] * a[i + 1][j];
        }
    }
    Gauss();
    for (int i = 1; i <= n; ++ i) printf("%.3lf ", d[i]);
    return 0;
}