P8141 [ICPC2020 WF] What’s Our Vector, Victor?

P8141 [ICPC2020 WF] What’s Our Vector, Victor?

给定 \(n\) 个 \(d\) 维向量 \(v_i\) 与 \(e_i\)，求 \(d\) 维向量 \(x\) 满足 \(\|x - v_i\| = e_i\)。保证有解。\(n, d \leq 500\)。

根据题意，可以列出以下方程。

\[\begin{cases} (x_1 - v_{1, 1}) ^ 2 + (x_2 - v_{1, 2}) ^ 2 + \cdots = e_1 ^ 2 \\ (x_1 - v_{2, 1}) ^ 2 + (x_2 - v_{2, 2}) ^ 2 + \cdots = e_2 ^ 2 \\ \cdots \end{cases} \]

将平方项展开后，出现了 \(x_i ^ 2\) 项。众所周知高斯消元只能处理线性方程组，所以考虑常用套路 差分消去平方项。处理后，我们得到了 \(n - 1\) 个 \(d\) 元线性方程组，和最开始的一个 \(d\) 元二次方程。形如

\[\begin{cases} (x_1 - v_{1, 1}) ^ 2 + (x_2 - v_{1, 2}) ^ 2 + \cdots = e_1 ^ 2 \\ 2(v_{1, 1} - v_{2, 1}) x_1 + 2(v_{1, 2} - v_{2, 2}) x_2 + \cdots = e_2 ^ 2 - e_1 ^ 2 + v_{1, 1} ^ 2 - v_{2, 1} ^ 2 + v_{1, 2} ^ 2 - v_{2, 1} ^ 2 + \cdots \\ 2(v_{2, 1} - v_{3, 1}) x_1 + 2(v_{2, 2} - v_{3, 2}) x_2 + \cdots = e_3 ^ 2 - e_2 ^ 2 + v_{2, 1} ^ 2 - v_{3, 1} ^ 2 + v_{2, 2} ^ 2 - v_{3, 1} ^ 2 + \cdots \\ \cdots \end{cases} \]

用向量和矩阵描述 \(x\) 前的系数，则有

\[\begin{cases} \| \pmb x - \pmb v_1 \| ^ 2 = e_1 ^ 2 \\ \pmb {Ax} = \pmb b \end{cases} \]

其中 \(\pmb A\) 为 \((n - 1) \times d\) 矩阵，\(\pmb x\) 为 \(d\) 维向量即所求答案，\(\pmb b\) 为 \(n - 1\) 维向量。

首先考虑化简二次式。令 \(\pmb {x'} = \pmb x - \pmb v_1\)。则方程可写为

\[\begin{cases} \| \pmb {x'} \| ^ 2 = e_1 ^ 2 \\ \pmb {A} \pmb x' = \pmb b - \pmb {Av}_1 \end{cases} \]

接下来我们用 \(\pmb {x'}\) 代替 \(\pmb x\)，求解后再加上 \(\pmb v_1\)。

高斯消元求解线性方程组 \(\pmb {A} \pmb x = \pmb b - \pmb {Av}_1\)，得到线性方程组的解集的 参数向量形式。即 \(\pmb x = \pmb x_{\rm base} + \sum\limits_{i = 1} ^ r c_i\pmb v_i\)，其中 \(r\) 是 自由元 个数。

• 众所周知，对于线性方程组的解集，任何 主元 都可以写作自由元的 线性组合 形式。考虑令所有自由元等于 \(0\)，我们得到了 \(\pmb x_{\rm base}\)。将所有自由元移到方程右侧，我们得到了所有主元关于所有自由元的线性组合的权，稍做变形即可得到对应向量。

  例如一个简化行阶梯型增广矩阵对应的线性方程组

  \[\begin{cases} x_1 + 6x_2 + 3x_4 = 0 \\ x_3 - 4x_4 = 5 \\ x_5 = 7 \end{cases} \]

  其中 \(x_1, x_3, x_5\) 是主元，\(x_2, x_4\) 是自由元。因此，其解集的参数表示为

  \[\begin{cases} x_1 = -6x_2 - 3x_4 \\ x_2 = x_2 \\ x_3 = 5 + 4x_4 \\ x_4 = x_4 \\ x_5 = 7 \end{cases} \]

  即 \(\pmb x = \begin{bmatrix} -6x_2 - 3x_4 \\ x_2 \\ 5 + 4x_4 \\ x_4 \\ 7 \end{bmatrix}\)。将常数，\(x_2\) 和 \(x_4\) 分离，得到 \(\pmb x = \begin{bmatrix} 0 \\ 0 \\ 5 \\ 0 \\ 7 \end{bmatrix} + x_2\begin{bmatrix} -6 \\ 1 \\ 0 \\ 0 \\ 0 \end{bmatrix} + x_4\begin{bmatrix} -3 \\ 0 \\ 4 \\ 1 \\ 0 \end{bmatrix}\)。即方程组的解可写为一个常向量 \(\pmb x_{\rm base}\) 加上若干向量的线性组合。线性组合的权对应了自由变量的取值。

为了让 \(\pmb x\) 满足 \(\| \pmb x \| = e_1\)，通常采用的手段是先求出 \(\pmb x\) 的 最小范数解 \(\pmb x_{\min}\)，再通过调整法求解。显然，如果 \(\| \pmb x_{\min}\| > e_1\) 则无解，与题意矛盾，因此 \(\| \pmb x_{\min} \| \leq e_1\)。而通过加上足够大的 \(\pmb v_i\)（任意一个均可），我们总能让得到的 \(\pmb x\) 的模长趋于无穷大。又因为 \(\| \pmb x \|\) 即 \(\| \pmb x_{\min} + k\pmb v_i \|\) 的变化随着 \(\pmb v_i\) 的权 \(k\) 的变化 连续，因此调整法正确。

• 关于最小范数解 \(\pmb x_{\min}\)，它与所有 \(\pmb v_i\) 垂直。若不与某个 \(\pmb v_i\) 垂直，令 \(\pmb x = \pmb x_{\min} + k\pmb v_i\) 且 \(\pmb x \perp \pmb v_i\)，通过点积的性质和数值分析容易证明 \(\| \pmb x \| < \| \pmb x_{\min} \|\)。因此可以列出方程 \(\pmb x_{\min} \cdot \pmb v_i = 0\ (1\leq i \leq r)\)，即 \(\pmb v_i ^ T \pmb x_{\min} = 0\)。若将 \(\pmb v_i\) 看成矩阵 \(\pmb V = \begin{bmatrix} \pmb v_1 & \pmb v_2 & \cdots & \pmb v_r \end{bmatrix}\)，则 \(\pmb V ^ T \pmb x_{\min} = \pmb 0\)。又因为 \(\pmb x = \pmb x_{\rm base} + \sum\limits_{i = 1} ^ r c_i \pmb v_i = \pmb x_{\rm base} + \pmb {Vc}\)，其中 \(\pmb c\) 是 \(r\times 1\) 的列向量，因此 \(\pmb V ^ T (\pmb x_{\rm base} + \pmb {Vc}) = \pmb 0\)。据此可高斯消元求得 \(\pmb c\)，继而求出对应的 \(\pmb x_{\min}\)。

由于 \(\pmb x_{\min}\) 垂直于所有 \(\pmb v_i\)，所以调整时直接使用勾股定理，随便选出一个向量 \(\pmb v_1\)，然后令 \(\pmb x = \pmb x_{\min} + \sqrt {e_1 ^ 2 - \| \pmb x_{\min} \| ^ 2} \times \dfrac {\pmb v_1} {\| \pmb v_1\|}\) 即可。时间复杂度 \(\mathcal{O}(n ^ 3)\)。

#include 
using namespace std;

const int N = 500 + 5;
const double eps = 1e-9;
int n, d, fr, ma, fre[N], mai[N], rev[N];
double v[N][N], A[N][N], dis[N], ans[N];
double b[N], xm[N], co[N][N];

int find(int m, double *A) {
	for(int i = 1; i <= m; i++) if(fabs(A[i]) > eps) return i;
	return -1;
}

template 
	void Gauss(int n, int m, T A) {
		int r = 1;
		for(int c = 1; c < m && r <= n; c++) {
			int p = r;
			for(int i = r + 1; i <= n; i++) if(fabs(A[i][c]) > fabs(A[p][c])) p = i;
			swap(A[r], A[p]);
			if(fabs(A[r][c]) < eps) continue;
			for(int i = m; i >= c; i--) A[r][i] /= A[r][c];
			for(int i = r + 1; i <= n; i++) for(int j = m; j >= c; j--) A[i][j] -= A[r][j] * A[i][c];
			r++;
		}
		
		for(int i = n - 1; i > 0; i--) {
			for(int j = i + 1; j <= n; j++) {
				int main = find(m - 1, A[j]);
				if(main == -1) continue;
				if(fabs(A[i][main]) > eps) for(int k = m; k >= main; k--) A[i][k] -= A[i][main] * A[j][k];
			}
			int main = find(m - 1, A[i]);
			if(main == -1) continue;
			for(int j = m; j >= main; j--) A[i][j] /= A[i][main];
		}
	}

void init() {
	for(int i = 1; i < n; i++) {
		int main = find(d, A[i]);
		if(main != -1) rev[main] = -1, mai[++ma] = main;
	}
	for(int i = 1; i <= d; i++) if(!rev[i]) fre[++fr] = i, co[rev[i] = fr][i] = 1;
}

void solve() {
	for(int i = 1; i < n; i++) {
		int main = find(d, A[i]);
		b[main] = A[i][d + 1];
		for(int j = main + 1; j <= d; j++) if(~rev[j]) co[rev[j]][main] = -A[i][j];
	}

	memset(A, 0, sizeof(A));
	for(int i = 1; i <= fr; i++) for(int j = 1; j <= d; j++) A[i][fr + 1] -= co[i][j] * b[j];
	for(int i = 1; i <= fr; i++) for(int j = 1; j <= d; j++) for(int k = 1; k <= fr; k++) A[i][k] += co[i][j] * co[k][j];
	Gauss(fr, fr + 1, A);
	for(int i = 1; i <= fr; i++) for(int j = 1; j <= d; j++) xm[j] += A[i][fr + 1] * co[i][j];
	for(int i = 1; i <= d; i++) xm[i] += b[i];
	
	double normx = 0, normc = 0;
	for(int i = 1; i <= d; i++) normx += xm[i] * xm[i], normc += co[1][i] * co[1][i];
	normc = sqrt(normc);
	for(int i = 1; i <= d; i++) ans[i] = xm[i] + co[1][i] / normc * sqrt(dis[1] * dis[1] - normx);
	
	
}

int main() {
	cin >> d >> n;
	for(int i = 1; i <= n; i++) {
		for(int j = 1; j <= d; j++) scanf("%lf", &v[i][j]);
		scanf("%lf", &dis[i]);
	}
	for(int i = 1; i < n; i++) {
		A[i][d + 1] = dis[i + 1] * dis[i + 1] - dis[i] * dis[i];
		for(int j = 1; j <= d; j++) {
			A[i][j] = 2 * (v[i][j] - v[i + 1][j]);
			A[i][d + 1] += v[i][j] * v[i][j] - v[i + 1][j] * v[i + 1][j];
		}
	}
	for(int i = 1; i < n; i++) for(int j = 1; j <= d; j++) A[i][d + 1] -= A[i][j] * v[1][j];
	
	Gauss(n - 1, d + 1, A), init();
	if(fr) solve();
	else for(int i = 1; i <= d; i++) ans[i] = A[i][d + 1];
	for(int i = 1; i <= d; i++) printf("%.17lf\n", ans[i] + v[1][i]);
	
	return 0;
}

/*
stupid mistakes:
	forget to check if main is -1
	注意 co[i][j] 的意义，表示 v_{i_{j, 1}}，即第 i 个自由元改变时，x_j 的变化系数 
	勾股定理忘记开根 =.=
	v1 正负搞错 
	求范数忘记开根 
*/