题目-球形空间产生器

在这里插入图片描述

问题分析

假设给定 n n n维, 设球心坐标 ( x 0 , x 1 , . . . , x n ) (x_0, x_1, ..., x_n) (x0,x1,...,xn), 假设给定点 1 1 1的坐标是 ( a 0 , a 1 , . . . , a n ) (a_0, a_1, ..., a_n) (a0,a1,...,an), 点 2 2 2的坐标是 ( b 0 , b 1 , . . . , b n ) (b_0, b_1, ..., b_n) (b0,b1,...,bn)

那么有
( a 0 − x 0 ) 2 + ( a 1 − x 1 ) 2 + . . . + ( a n − x n ) 2 = ( b 0 − x 0 ) 2 + ( b 1 − x 1 ) 2 + . . . + ( b n − x n ) 2 (a_0 - x_0) ^ 2 + (a_1 - x_1) ^ 2 + ... + (a_n - x_n) ^ 2 = (b_0 - x_0) ^ 2 + (b_1 - x_1) ^ 2 + ... + (b_n - x_n) ^ 2 (a0x0)2+(a1x1)2+...+(anxn)2=(b0x0)2+(b1x1)2+...+(bnxn)2

化简得到

a 0 2 − 2 a 0 x 0 + a 1 2 − 2 a 1 x 1 + . . . + a n 2 − 2 a n x n = b 0 2 − 2 b 0 x 0 + b 1 2 − 2 b 1 x 1 + . . . . + b n 2 − 2 b n x n a_0 ^ 2 - 2a_0x_0 + a_1 ^ 2 - 2a_1x_1 + ... + a_n ^ 2 - 2a_nx_n = b_0 ^ 2 - 2b_0x_0 + b_1 ^ 2 - 2b_1x_1 + .... + b_n ^ 2 - 2b_nx_n a022a0x0+a122a1x1+...+an22anxn=b022b0x0+b122b1x1+....+bn22bnxn

继续化简

2 ( a 0 − b 0 ) ⋅ x 0 + 2 ( a 1 − b 1 ) ⋅ x 1 + . . . . + 2 ( a n − b n ) ⋅ x n = a 0 2 + a 1 2 + . . . + a n 2 − b 0 2 − b 1 2 . . . − b n 2 2(a_0 - b_0) \cdot x_0 + 2(a_1 - b_1) \cdot x_1 + .... + 2(a_n - b_n) \cdot x_n = a_0 ^ 2 + a_1 ^ 2 + ... + a_n ^ 2 - b_0 ^ 2 - b_1 ^ 2 ... - b_n ^ 2 2(a0b0)x0+2(a1b1)x1+....+2(anbn)xn=a02+a12+...+an2b02b12...bn2

也就是
( a 0 − b 0 ) ⋅ x 0 + ( a 1 − b 1 ) ⋅ x 1 + . . . . + ( a n − b n ) ⋅ x n = a 0 2 + a 1 2 + . . . + a n 2 − b 0 2 − b 1 2 . . . − b n 2 2 (a_0 - b_0) \cdot x_0 + (a_1 - b_1) \cdot x_1 + .... + (a_n - b_n) \cdot x_n = \frac{a_0 ^ 2 + a_1 ^ 2 + ... + a_n ^ 2 - b_0 ^ 2 - b_1 ^ 2 ... - b_n ^ 2}{2} (a0b0)x0+(a1b1)x1+....+(anbn)xn=2a02+a12+...+an2b02b12...bn2

两个坐标能得到一个方程, 给定 n + 1 n + 1 n+1个坐标, 那么能得到 n n n个方程, n n n个方程, n n n个未知数, 用高斯消元算法求解

算法步骤

  • 构建线性方程组
  • 用高斯消元算法求解

算法时间复杂度 O ( n 3 ) O(n ^ 3) O(n3)

代码实现

#include <bits/stdc++.h>

using namespace std;

const int N = 15;
const double EPS = 1e-10;

int n;
double m[N][N];

void gauss() {
   
    int r, c;
    for (r = c = 0; c < n; ++c) {
   
        int idx = r;
        for (int i = r + 1; i < n; ++i) {
   
            if (fabs(m[i][c]) > fabs(m[idx][c])) idx = i;
        }

        if (fabs(m[idx][c]) < EPS) continue;

        for (int i = c; i <= n; ++i) swap(m[r][i], m[idx][i]);
        for (int i = n; i >= c; --i) m[r][i] /= m[r][c];
        for (int i = r + 1; i < n; ++i) {
   
            for (int j = n; j >= c; --j) {
   
                m[i][j] -= m[r][j] * m[i][c];
            }
        }
        r++;
    }

    for (int i = n - 1; i >= 0; --i) {
   
        for (int j = i + 1; j < n; ++j) {
   
            m[i][n] -= m[i][j] * m[j][n];
        }
    }

    for (int i = 0; i < n; ++i) printf("%.3lf ", m[i][n]);
}

int main() {
   
    ios::sync_with_stdio(false);
    cin.tie(0);

    cin >> n;

    double tmp[N] = {
   0};
    for (int i = 0; i < n; ++i) cin >> tmp[i];

    for (int i = 0; i < n; ++i) {
   
        for (int j = 0; j < n; ++j) {
   
            double val;
            cin >> val;

            m[i][j] = val - tmp[j];
            m[i][n] += (val * val - tmp[j] * tmp[j]) / 2;
            tmp[j] = val;
        }
    }

    gauss();

    return 0;
}