P8141 [ICPC 2020 WF] What’s Our Vector, Victor? 题解（ChatGPT）

• https://www.zhihu.com/question/581556221/answer/97377273573
• [ICPC 2020 WF] What’s Our Vector, Victor?

题目理解

题目要求在 \(d\) 维空间中确定一个未知点 \(x\)（你的位置），已知空间中 \(n\) 个蜗牛的坐标 \(p_1, p_2, \dots, p_n \in \mathbb{R}^d\) 以及你与每只蜗牛之间的欧几里得距离 \(r_1, r_2, \dots, r_n\)（数据保证一致，即存在一个 \(x\) 同时满足所有距离方程）。数学上，对于所有 \(i=1,\dots,n\) 有

\[\|x - p_i\| = r_i. \]

由于数据均由真实位置生成，所以这些方程是一致的。

---

解题思路

1. 消去二次项构造线性方程组

对每个 \(i\)，我们有

\[\|x-p_i\|^2 = r_i^2. \]

选取第一个蜗牛作为参考，将 \(i\ge2\) 的方程与 \(i=1\) 的方程相减：

\[\|x-p_i\|^2 - \|x-p_1\|^2 = r_i^2 - r_1^2. \]

展开后，注意到 \(x^T x\) 相互抵消，得到

\[-2p_i^T x + \|p_i\|^2 + 2p_1^T x - \|p_1\|^2 = r_i^2 - r_1^2, \]

整理后写为

\[(p_i - p_1)^T x = \frac{r_1^2 - r_i^2 + \|p_i\|^2 - \|p_1\|^2}{2}. \]

记

\[A_i = p_i - p_1,\quad b_i = \frac{r_1^2 - r_i^2 + \|p_i\|^2 - \|p_1\|^2}{2},\quad (i=2,\dots,n). \]

这样就得到了一个 $ (n-1) \times d$ 的线性方程组

\[A x = b. \]

2. 求解线性系统

• 满秩情况：若方程组满秩（通常当 \(n\ge d+1\) 且点的分布足够“通用”时成立），则线性系统有唯一解。用高斯消元求得特解 \(x_0\)（对自由变量直接置零），理论上应满足所有差分方程，并且利用第一个方程（\(\|x-p_1\|=r_1\)）也能成立。

• 欠定情况或数值误差：当 \(n-1 < d\) 或矩阵 \(A\) 的秩小于 \(d\) 时，系统欠定，得到的特解 \(x_0\)可能不满足第一个方程，即可能有

  \[\|x_0-p_1\| \neq r_1. \]

  但由于数据一致，实际存在解。注意到一般解可写成

  \[x = x_0 + v,\quad \text{其中 } v \in \mathrm{Null}(A). \]

  为了修正 \(x_0\) 使得 \(\|x-p_1\| = r_1\)，我们作如下处理：

  • 记 \(w_0 = x_0 - p_1\)。将 \(w_0\) 分解为行空间和零空间两部分。利用对 \(A\) 的行空间做 Gram–Schmidt 正交化，计算 \(w_0\) 在行空间上的投影，然后零空间分量即为
    \[v = w_0 - \mathrm{Proj}_{\mathrm{Row}(A)}(w_0). \]

  • 令修正解为
    \[x = x_0 + \alpha v, \]
    并要求
    \[\|x-p_1\| = \|w_0 + \alpha v\| = r_1. \]
    展开后得到一元二次方程（关于 \(\alpha\)），求解后选择合适的 \(\alpha\) 使得调整最小，从而获得满足所有条件的解。

• 满秩但由于数值原因误差较大时：若系统满秩但得到的 \(x_0\) 与参考距离 \(r_1\) 仍有较大误差，可以沿着 \(p_1 \to x_0\) 的方向进行缩放修正，使得最终 \(\|x-p_1\|=r_1\)。

• 特殊情况 \(n=1\)：只有一个距离信息时，任取 \(p_1 + (r_1,0,\dots,0)\) 即可。

3. 数值精度

由于题目中数据以 17 位精度给出，且 \(d, n\) 的范围最多为 500，为了保证误差不超过 \(10^{-5}\)，所有计算均采用 long double 类型，并且在高斯消元、正交化等步骤中设置了较小的容差。

---

代码时间复杂度分析

代码主要分为以下几个部分：

1. 高斯消元求解线性系统

   • 对一个 \(m \times d\) 的矩阵（其中 \(m = n-1\)）进行消元，复杂度大致为 \(O(m \cdot d^2)\)。
   • 最坏情况下 \(n, d \le 500\)，故该部分的复杂度为 \(O(500 \times 500^2) = O(125\,000\,000)\) 次基本运算。

2. Gram–Schmidt 正交化

   • 用于计算 \(A\) 的行空间正交基。设矩阵行数为 \(m\)，每个向量维数为 \(d\)，此部分的复杂度为 \(O(m \cdot d^2)\)。
   • 最坏情况同样为 \(O(125\,000\,000)\) 次基本运算。

3. 其他向量加减、点积和范数计算

   • 每次操作均为 \(O(d)\)，且总次数不会超过 \(O(n \cdot d)\) 或 \(O(m \cdot d)\)，在本题范围内开销可以忽略。

综上，总的时间复杂度为

\[O(m \cdot d^2) = O(n \cdot d^2) \]

在最坏情况下 \(n = 500,\, d = 500\) 时约为 \(125\)M 次基本运算，加上常数因子（以及 long double 可能的额外开销），但由于数据规模较小且常数较低，该算法在 C++ 编译器优化（如 -O2）下能够在 1000ms 时限内顺利通过。

---

总结

本题利用“多边定位”（multilateration）的思想，通过对距离平方方程做差分消去二次项，构造出一个线性方程组，从而求得候选解 \(x_0\)。当系统欠定或数值误差较大时，通过将 \(x_0\) 在零空间上做适当调整，使得最终解满足所有距离条件。整个过程使用高精度 long double 保证结果精度，并且时间复杂度为 \(O(n \cdot d^2)\) ，能够满足题目数据范围要求。

代码：

#include <bits/stdc++.h>
using namespace std;// 使用 long double 以提高计算精度
const long double EPS = 1e-18L;
const long double TOL = 1e-7L;// 计算向量 v 的欧式范数
long double normVec(const vector<long double>& v) {long double s = 0.0L;for (long double x : v) s += x * x;return sqrt(s);
}// 计算两个向量的点积
long double dotVec(const vector<long double>& a, const vector<long double>& b) {long double s = 0.0L;int n = a.size();for (int i = 0; i < n; i++) s += a[i] * b[i];return s;
}// 高斯消元：求解 A*x = b
// A 为 m×d 矩阵，b 为 m 维向量。
// 在欠定系统中，令自由元为 0。
// 同时输出矩阵的秩以及主元所在列（pivot_col）。
vector<long double> gaussElim(vector<vector<long double>> A, vector<long double> b, int d, int m, int &rank, vector<int> &pivot_col) {rank = 0;int row = 0, col = 0;pivot_col.clear();while(row < m && col < d) {// 在 col 列中从 row 到 m-1 寻找主元int pivot = row;for (int i = row; i < m; i++){if (fabsl(A[i][col]) > fabsl(A[pivot][col]))pivot = i;}if(fabsl(A[pivot][col]) < EPS){col++;continue;}if(pivot != row){swap(A[row], A[pivot]);swap(b[row], b[pivot]);}pivot_col.push_back(col);long double pivVal = A[row][col];for (int j = col; j < d; j++){A[row][j] /= pivVal;}b[row] /= pivVal;for (int i = 0; i < m; i++){if(i != row){long double factor = A[i][col];if(fabsl(factor) > EPS){for (int j = col; j < d; j++){A[i][j] -= factor * A[row][j];}b[i] -= factor * b[row];}}}row++;col++;}rank = row;vector<long double> x(d, 0.0L);for (int i = 0; i < rank; i++){int c = pivot_col[i];x[c] = b[i];}return x;
}// 用 Gram–Schmidt 正交化法计算矩阵 A（每一行视为 R^d 中的一个向量）的行空间正交基
vector<vector<long double>> computeOrthonormalBasis(const vector<vector<long double>> &A) {int m = A.size();int d = (m ? A[0].size() : 0);vector<vector<long double>> basis;for (int i = 0; i < m; i++){vector<long double> v = A[i];for(auto &q : basis) {long double proj = dotVec(v,q);for (int j = 0; j < d; j++){v[j] -= proj * q[j];}}long double normv = normVec(v);if(normv > EPS){for (int j = 0; j < d; j++){v[j] /= normv;}basis.push_back(v);}}return basis;
}// 主函数
int main(){ios::sync_with_stdio(false);cin.tie(nullptr);int d, n;cin >> d >> n;vector<vector<long double>> p(n, vector<long double>(d, 0.0L));vector<long double> r(n, 0.0L);for (int i = 0; i < n; i++){for (int j = 0; j < d; j++){double tmp;cin >> tmp;p[i][j] = tmp;}double tmp;cin >> tmp;r[i] = tmp;}// 若只有一个蜗牛，则随便选一个在 p[0] 距离 r[0] 的点if(n == 1){vector<long double> ans = p[0];ans[0] += r[0];cout << fixed << setprecision(12);for (int i = 0; i < d; i++){cout << (long double)ans[i] << ( i==d-1 ? "\n" : " ");}return 0;}// 以第一个蜗牛的位置 p[0] 为参考vector<long double> p1 = p[0];long double r1 = r[0];// 对 i=2..n 建立差分方程： (p[i]-p1)^T x = (r1^2 - r[i]^2 + ||p[i]||^2 - ||p1||^2)/2int m = n - 1;vector<vector<long double>> A_orig(m, vector<long double>(d, 0.0L));for (int i = 1; i < n; i++){for (int j = 0; j < d; j++){A_orig[i-1][j] = p[i][j] - p1[j];}}// 为消元保存一份副本vector<vector<long double>> A = A_orig;vector<long double> b(m, 0.0L);long double normP1Sq = 0.0L;for (int j = 0; j < d; j++){normP1Sq += p1[j]*p1[j];}for (int i = 1; i < n; i++){long double normPiSq = 0.0L;for (int j = 0; j < d; j++){normPiSq += p[i][j]*p[i][j];}long double r1sq = r1*r1, ri_sq = r[i]*r[i];b[i-1] = (r1sq - ri_sq + normPiSq - normP1Sq) / 2.0L;}int rank;vector<int> pivot_col;vector<long double> x0 = gaussElim(A, b, d, m, rank, pivot_col);// 检查 x0 是否满足第一个蜗牛的距离条件： ||x0-p1|| = r1vector<long double> diff(d, 0.0L);for (int j = 0; j < d; j++){diff[j] = x0[j] - p1[j];}long double norm_diff = normVec(diff);long double err = fabsl(norm_diff - r1);if(err < TOL){cout << fixed << setprecision(12);for (int j = 0; j < d; j++){cout << (long double)x0[j] << (j==d-1 ? "\n" : " ");}return 0;}// 如果系统欠定（rank < d），则利用零空间进行修正：// 记 w0 = x0-p1，将其在 A_orig 行空间上的正交投影记为 proj，// 则零空间成分为 v = w0 - proj. 令修正后解为 x_new = x0 + α*v，// 要求 ||(x0 + α*v) - p1|| = r1，从而解出 α。if(rank < d) {vector<vector<long double>> Q = computeOrthonormalBasis(A_orig); // 行空间正交基vector<long double> w0 = diff; // = x0-p1vector<long double> proj(d, 0.0L);for(auto &q : Q) {long double dp = dotVec(w0, q);for (int j = 0; j < d; j++){proj[j] += dp * q[j];}}// v 为 w0 在零空间中的分量vector<long double> v(d, 0.0L);for (int j = 0; j < d; j++){v[j] = w0[j] - proj[j];}long double norm_v = normVec(v);if(norm_v < EPS){// 极少数情况下无可修正分量，输出 x0cout << fixed << setprecision(12);for (int j = 0; j < d; j++){cout << (long double)x0[j] << (j==d-1 ? "\n" : " ");}return 0;}// 求解 α，使得 ||w0 + α*v|| = r1// 即: ||w0||^2 + 2α (w0·v) + α^2||v||^2 = r1^2.long double a_coef = norm_v * norm_v;long double b_coef = 2.0L * dotVec(w0, v);long double c_coef = norm_diff * norm_diff - r1 * r1;long double disc = b_coef * b_coef - 4.0L * a_coef * c_coef;if(disc < 0) disc = 0;long double sqrt_disc = sqrt(disc);long double alpha1 = (-b_coef + sqrt_disc) / (2.0L * a_coef);long double alpha2 = (-b_coef - sqrt_disc) / (2.0L * a_coef);long double alpha = (fabsl(alpha1) < fabsl(alpha2) ? alpha1 : alpha2);vector<long double> x_new(d, 0.0L);for (int j = 0; j < d; j++){x_new[j] = x0[j] + alpha * v[j];}// 检查修正后是否满足要求vector<long double> diff_new(d, 0.0L);for (int j = 0; j < d; j++){diff_new[j] = x_new[j] - p1[j];}long double norm_new = normVec(diff_new);if(fabsl(norm_new - r1) < TOL){cout << fixed << setprecision(12);for (int j = 0; j < d; j++){cout << (long double)x_new[j] << (j==d-1 ? "\n" : " ");}return 0;} else {// 若仍有微小误差，直接输出修正解cout << fixed << setprecision(12);for (int j = 0; j < d; j++){cout << (long double)x_new[j] << (j==d-1 ? "\n" : " ");}return 0;}}// 若系统满秩（rank == d）而误差仍较大，则说明由于数值问题造成误差（理论上唯一解应同时满足所有条件）。// 此时我们用一个一维牛顿法沿 p1->x0 的射线修正（注意：这种调整会破坏差分方程，但数据本身一定一致）。long double t = r1 / norm_diff;vector<long double> x_new(d, 0.0L);for (int j = 0; j < d; j++){x_new[j] = p1[j] + t*(x0[j]-p1[j]);}cout << fixed << setprecision(12);for (int j = 0; j < d; j++){cout << (long double)x_new[j] << (j==d-1 ? "\n" : " ");}return 0;
}