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\) 使得调整最小，从而获得满足所有条件的解。
- 记 \(w_0 = x_0 - p_1\)。将 \(w_0\) 分解为行空间和零空间两部分。利用对 \(A\) 的行空间做 Gram–Schmidt 正交化，计算 \(w_0\) 在行空间上的投影，然后零空间分量即为
-
满秩但由于数值原因误差较大时：若系统满秩但得到的 \(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 类型，并且在高斯消元、正交化等步骤中设置了较小的容差。
代码时间复杂度分析
代码主要分为以下几个部分：
-
高斯消元求解线性系统
- 对一个 \(m \times d\) 的矩阵（其中 \(m = n-1\)）进行消元，复杂度大致为 \(O(m \cdot d^2)\)。
- 最坏情况下 \(n, d \le 500\)，故该部分的复杂度为 \(O(500 \times 500^2) = O(125\,000\,000)\) 次基本运算。
-
Gram–Schmidt 正交化
- 用于计算 \(A\) 的行空间正交基。设矩阵行数为 \(m\)，每个向量维数为 \(d\)，此部分的复杂度为 \(O(m \cdot d^2)\)。
- 最坏情况同样为 \(O(125\,000\,000)\) 次基本运算。
-
其他向量加减、点积和范数计算
- 每次操作均为 \(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)/2
int 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|| = r1
vector<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-p1
vector<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){
// 极少数情况下无可修正分量，输出 x0
cout << 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;
}
浙公网安备 33010602011771号