LOJ 6803 - 「ICPC World Finals 2020」胜利者，我们的向量是什么？（高斯消元）

LOJ 题面传送门

提供一种略有不同的做法。思路来自 ljc1301。

首先假设我们要求的答案为 \(x_1,x_2,\cdots,x_d\)，给定的坐标为 \((a_{i,1},a_{i,2},\cdots,a_{i,d})\)，待求坐标离这 \(n\) 个坐标的距离分别为 \(D_1,D_2,\cdots,D_n\)，那么很显然我们可以列出 \(n\) 个方程组

\[\begin{cases} (a_{1,1}-x_1)^2+(a_{1,2}-x_2)^2+\cdots+(a_{1,d}-x_d)^2=D_1^2\\ (a_{2,1}-x_1)^2+(a_{2,2}-x_2)^2+\cdots+(a_{2,d}-x_d)^2=D_2^2\\ \cdots\\ (a_{n,1}-x_1)^2+(a_{n,2}-x_2)^2+\cdots+(a_{n,d}-x_d)^2=D_n^2 \end{cases} \]

方程中待二次，不好直接计算，因此考虑做差，设 \(b_{i,j}=2(a_{i,j}-a_{i+1,j})\)，\(C_i=D_{i+1}^2-D_i^2+\sum\limits_{t=1}^da_{i,t}^2-a_{i+1,t}^2\)。那么一组 \(\{x_d\}\) 合法的必要条件就是

\[\begin{cases} b_{1,1}x_1+b_{1,2}x_2+\cdots+b_{1,d}x_d=C_1\\ b_{2,1}x_1+b_{2,2}x_2+\cdots+b_{2,d}x_d=C_2\\ \cdots\\ b_{n-1,1}x_1+b_{n-1,2}x_2+\cdots+b_{n-1,d}x_d=C_{n-1} \end{cases} \]

根据线性代数那套理论，设 \(m\) 表示该方程组的主元个数，\(y_1,y_2,\cdots,y_m\) 表示方程组的 \(m\) 个主元，那么所有的 \(x_i\) 均可表示为 \(y_1,y_2,\cdots,y_m\) 的线性组合，这样我们只要解出任意一组 \(y_i\) 即可知道对应的 \(x_i\)。

不过别忘了，\(\{x_d\}\) 是上述方程组的解只是这组 \(\{x_d\}\) 合法的必要条件，要使得 \(\{x_d\}\) 真正合法，还要有

\[(a_{1,1}-x_1)^2+(a_{1,2}-x_2)^2+\cdots+(a_{1,d}-x_d)^2=D_1^2 \]

我们将所有 \(x_i\) 替换为 \(y_i\) 的线性组合，那么等式就可以写成

\[\begin{aligned} Q_{1,1}y_1^2+Q_{1,2}y_1y_2+Q_{1,3}y_1y_3\cdots+Q_{1,m}y_1y_m+Q_{1,m+1}y_1+&\\ Q_{2,2}y_2^2+Q_{2,3}y_2y_3\cdots+Q_{2,m}y_2y_m+Q_{2,m+1}y_2+&\\ \cdots&\\ Q_{m,m}y_m^2+Q_{m,m+1}y_m+&\\ Q_{m+1,m+1}&\\ =0& \end{aligned} \]

然后就是比较关键的一步了，考虑主元配方，即将含 \(y_1\) 的项（\(Q_{1,1}y_1^2,Q_{1,i}y_1y_i,Q_{1,m+1}y_1\)）的项配方，这样含有 \(y_1\) 的部分就可以写成 \(\square(y_1+\square y_2+\square y_3+\cdots+\square y_m+\square)^2\) 的形式，其中 \(\square\) 均代表一个我们可以求得的常数，配完了 \(y_1\)，接下来以 \(y_2\) 为主元，将剩余部分与 \(y_2\) 有关的交叉项都配成 \(\square(y_2+\square y_3+\square y_4+\cdots+\square y_m+\square)^2\)​ 的形式，依此类推，这样最终的式子可以写成

\[\begin{aligned} \square(y_1+\square y_2+\square y_3+\square y_4+\cdots+\square y_{m-1}+\square y_m+\square)^2+&\\ \square(y_2+\square y_3+\square y_4+\cdots+\square y_{m-1}+\square y_m+\square)^2+&\\ \cdots+&\\ \square(y_{m-1}+\square y_m+\square)^2+&\\ \square(y_m+\square)^2=&\\ V \end{aligned} \]

由于题目保证有解，所以必然 \(V\ge 0\)，此时我们就令最后一项结果为 \(V\)，解出 \(y_m\)，再根据前面所有的括号都为 \(0\) 的条件，顺着推回去即可。

时间复杂度 \(n^3\)。

const int MAXN = 500;
const ld EPS = 1e-11;
int n, d, t[MAXN + 5], m, id[MAXN + 5];
ld x[MAXN + 5][MAXN + 5], D[MAXN + 5], A[MAXN + 5][MAXN + 5], B[MAXN + 5][MAXN + 5];
int main() {
//	freopen("snail.in", "r", stdin);
//	freopen("snail.out", "w", stdout);
	scanf("%d%d", &d, &n); srand(time(0));
	for (int i = 1; i <= n; i++) {
		for (int j = 1; j <= d; j++) scanf("%Lf", &x[i][j]);
		scanf("%Lf", &D[i]);
	}
	for (int i = 1; i < n; i++) {
		for (int j = 1; j <= d; j++) A[i][j] = 2.0 * (x[i][j] - x[i + 1][j]);
		A[i][d + 1] = D[i + 1] * D[i + 1] - D[i] * D[i];
		for (int j = 1; j <= d; j++) {A[i][d + 1] += x[i][j] * x[i][j]; A[i][d + 1] -= x[i + 1][j] * x[i + 1][j];}
	}
	int cur = 1;
	for (int i = 1; i <= d; i++) {
		int pos = 0;
		for (int k = cur; k < n; k++) if (fabs(A[k][i]) > EPS) pos = k;
		if (pos) {
			for (int k = i; k <= d + 1; k++) swap(A[pos][k], A[cur][k]);
			for (int k = i + 1; k <= d + 1; k++) A[cur][k] /= A[cur][i];
			A[cur][i] = 1;
			for (int k = cur + 1; k < n; k++) {
				for (int l = i + 1; l <= d + 1; l++) A[k][l] -= A[cur][l] * A[k][i];
				A[k][i] = 0;
			}
			++cur;
		} else t[++m] = i, id[i] = m;
	}
	if (m == 0) {
		static ld res[MAXN + 5];
		for (int i = cur - 1; i; i--) {
			res[i] = A[i][d + 1];
			for (int j = i + 1; j <= d; j++) res[i] -= res[j] * A[i][j];
		}
		for (int i = 1; i <= d; i++) printf("%.12Lf%c", res[i], " \n"[i == d]);
		return 0;
	}
	for (int i = cur - 1; i; i--) {
		int pos = 0;
		for (int j = 1; j <= d; j++) if (fabs(A[i][j]) > 0) {
			pos = j; break;
		}
		for (int j = pos + 1; j <= d; j++) {
			if (id[j]) B[pos][id[j]] += -A[i][j];
			else {
				for (int k = 1; k <= m + 1; k++)
					B[pos][k] -= B[j][k] * A[i][j];
			}
		}
		B[pos][m + 1] += A[i][d + 1];
	}
	for (int i = 1; i <= d; i++) if (id[i]) B[i][id[i]] = 1;
	static ld coef[MAXN + 5][MAXN + 5];
	for (int i = 1; i <= d; i++) {
		for (int j = 1; j <= m; j++) for (int k = j + 1; k <= m; k++)
			coef[j][k] += 2.0 * B[i][j] * B[i][k];
		for (int j = 1; j <= m; j++) coef[j][j] += B[i][j] * B[i][j];
		for (int j = 1; j <= m; j++) coef[j][m + 1] += 2.0 * (B[i][m + 1] - x[1][i]) * B[i][j];
		coef[m + 1][m + 1] += (B[i][m + 1] - x[1][i]) * (B[i][m + 1] - x[1][i]);
	}
	coef[m + 1][m + 1] -= D[1] * D[1];
	static ld eqs[MAXN + 5][MAXN + 5], coef_before[MAXN + 5];
	for (int i = 1; i <= m; i++) {
		ld pre = coef[i][i];
		for (int j = i; j <= m + 1; j++) coef[i][j] /= pre;
		eqs[i][i] = 1; coef_before[i] = pre;
		for (int j = i + 1; j <= m + 1; j++) eqs[i][j] = coef[i][j] / 2.0;
		for (int j = i + 1; j <= m + 1; j++) for (int k = i + 1; k <= m + 1; k++)
			coef[min(j, k)][max(j, k)] -= pre * eqs[i][j] * eqs[i][k];
	}
	static ld res[MAXN + 5];
	ld V = -coef[m + 1][m + 1] / coef_before[m];
	V = sqrt(V); res[m] = (V - eqs[m][m + 1]) / eqs[m][m];
	res[m + 1] = 1;
	for (int i = m - 1; i; i--) {
		res[i] = 0;
		for (int j = i + 1; j <= m + 1; j++) res[i] -= res[j] * eqs[i][j];
		res[i] /= eqs[i][i];
	}
	for (int i = 1; i <= d; i++) {
		ld sum = 0;
		for (int j = 1; j <= m + 1; j++) sum += res[j] * B[i][j];
		printf("%.10Lf%c", sum, " \n"[i == d]);
	}
	return 0;
}