一、引言

如果现在你想制作一张二维地图，上面要标出西安的兵马俑、北京的故宫、上海的东方明珠，还有苏州的每一座园林、武汉的每一家热干面馆，这做得到吗？

有人说，做得到，顺便打开了高德地图或者百度地图。可是别忘了，地球是个球体，上面的这些地点其实都位于球面上。所以所有的二维平面地图都有一定程度的失真 —— 比如俄罗斯在地图上似乎有中国的三四倍大小，但实际面积却只有 1.7 倍。

好吧，那就做个三维地图就好了？就像地球仪那样。这样确实有道理，但想像这么一个场景：你想知道从北京的故宫沿直线走路去上海的东方明珠所需要的时间，于是你找到了它们在三维地图中分别对应的两个点 —— 下一步，是测量下这两个点的直线距离，然后除以比例尺，再除以走路的速度吗？

答案是否定的！

因为你走路只会在球面上走，所以你所谓的“沿直线走路”，并不是真的直线，而是球面上连接两点的弧线中最短的那条。所以球面上两点之间的“距离”，其本质乃是球面上一段弧的长度，而非在三维空间中直接沿直线相连的长度。要想从三维地图上知道你走路所需要的时间，你需要拿一把尺子，先把零刻度线对准北京故宫，再向着上海东方明珠的方向，缓缓沿球面转动尺子，直到某个刻度线碰到了上海的东方明珠；这时你把此时尺子的读数除以比例尺，再除以走路速度，这才是你想计算的时间。

那么问题来了 —— 给定一个球面，你能否找到一种办法，把球面上的每个点都唯一地映射到三维空间中，使得每个点原本在球面上的距离，恰好等于映射后在三维空间中的直线距离？如果三维不行，四维呢？五维呢？$n$ 维呢？

这正是本次 UER #10 的 B 题 重构元宇宙 想讨论的问题。该题大意是：给定 $n$ 个点两两之间的距离，你能把他们画在 $k$ 维空间中，且保持所有点两两之间的距离不变吗？还有个附加条件是，对两点之间的距离的查询次数不能超过 $n(k+1)$。为避免精度误差，原题的点的坐标是随机的。这里还有一个加强版，去掉了这个随机性的条件，供勇士们挑战。

所以，要想知道球面上的距离能不能被映射为某个欧式空间的距离，我们只需要先假设它们可以表示为 $k$ 维空间中的距离，然后带到这道 B 题的程序里，看看程序输出的结果是否真的满足假设。

实际生活中或者科学研究中，我们会碰到很多种“距离”，它们的具体定义取决于实际的应用，而非总是像欧式空间一样，把每一维坐标相减，平方加起来，再开根号。例如上面举例的球面上的距离，人们是为了地理上的用途而定义它，在计算中国到美国的飞机行程时也会使用它。如果有一个想从中国去往美国的旅行者，对欧式距离爱得深沉，一定要让自己旅行的时间与欧式距离成正比，那他只能强行通过钻入地下、穿过地心的方式去往美国了。

那么这道题到底该怎么做呢？

二、转换为向量内积

首先，距离其实不是个很好操作的东西。但是内积就不一样了。关于距离的问题，有时候可以转化为关于内积的问题，就像做最小圆覆盖时求三角形外接圆的方法一样。（这一步 hehezhou 已经讲得很好了，不过我还是再写一遍 hhh）

由于平移整个点集不会影响两两之间的距离，我们不妨设 $0$ 号点坐标为 $0$，并以此为原点建立一个 $k$ 维的直角坐标系。对于 $1$ 到 $n-1$ 号点，设 $i$ 号点第 $j$ 维的坐标为 $x_{i,j}$，那么已知 $p$ 号点和 $q$ 号点之间的距离为 $d(p, q)$，实际意思是： $$ \sum_{j=1}^{k} (x_{p,j} - x_{q,j})^2 = d^2(p, q). $$ 注意我们已经把 $0$ 号点作为原点了，所以如果我们知道 $p$ 号点和 $0$ 号点之间的距离为 $d(p, 0)$，那么： $$ \sum_{j=1}^{k} x_{p,j}^2 = d^2(p, 0). $$ 同样，对于 $q$ 号点我们也有： $$ \sum_{j=1}^{k} x_{q,j}^2 = d^2(q, 0). $$ 现在把第一个式子减去后两个式子，就可以得到：

\begin{align} \sum_{j=1}^{k}\left( (x_{p,j} - x_{q,j})^2 - x_{p,j}^2 - x_{q,j}^2 \right) &= d^2(p, q) - d^2(p, 0) - d^2(q, 0) \\ \sum_{j=1}^{k} 2 x_{p,j} \cdot x_{q,j} &= d^2(p, q) - d^2(p, 0) - d^2(q, 0) \\ \sum_{j=1}^{k} x_{p,j} \cdot x_{q,j} &= \frac{1}{2}\left(d^2(p, q) - d^2(p, 0) - d^2(q, 0) \right). \end{align}

令 $A_{p,q} = \frac{1}{2}\left(d^2(p, q) - d^2(p, 0) - d^2(q, 0) \right)$，那么原问题就变成了一个已知 $n - 1$ 个向量两两之间的内积，构造一组满足条件的向量坐标的问题了。 $$ \forall p, q \in \{1, \dots, n -1\}: \qquad \sum_{j=1}^{k} x_{p,j} \cdot x_{q,j} = A_{p,q}. \tag{F} $$ 大家可能在高中数学必修四里学过平面向量的内积，高维是类似的，都是将两个向量每个维度的坐标相乘再相加。上式左边就是 $x_p$ 和 $x_q$ 这两个向量的内积。

我们可以通过 $n-1$ 次查询得到每一个点到原点的距离 $d(\,\cdot\,, 0)$，然后每次想知道 $A_{p,q}$ 时，只要额外查询一下 $d(p, q)$ 就好了。

三、高斯消元？

hehezhou 的题解已经介绍了如何使用高斯消元解决本题，这里不再赘述。但是如果数据不是随机数据，那么高斯消元的误差可能会非常大！可以思考一下下面这个例子（也许我之后会再写一篇怎么卡高斯消元的博客，并且解释下为啥随机数据不会出问题233）：

11 10 121
0 0 0 0 0 0 0 0 0 0
4 1 0 0 0 0 0 0 0 0
-224 -57 1 0 0 0 0 0 0 0
-224 0 -57 1 0 0 0 0 0 0
-224 0 0 -57 1 0 0 0 0 0
-224 0 0 0 -57 1 0 0 0 0
-224 0 0 0 0 -57 1 0 0 0
-224 0 0 0 0 0 -57 1 0 0
-224 0 0 0 0 0 0 -57 1 0
-224 0 0 0 0 0 0 0 -57 1
-227 0 0 0 0 0 0 0 0 -57

四、楚列斯基分解

我想补充一个常数更小、精度更高的算法：楚列斯基分解（Cholesky decomposition）。这个算法想要解决的问题，就是求出一组满足 (F) 的一组向量（其实还会额外保证 $x_{p,j}$ 在 $j > p$ 的时候等于 $0$）。使用这个算法不仅能通过原 B 题，还能通过加强版。

楚列斯基是个法国数学家，全名 André-Louis Cholesky，出生于 1875 年。在从事大地测量学和制图学工作的时候，他想到了这个楚列斯基分解算法。但是后来第一次世界大战爆发，楚列斯基就从军去了，成为了一名法国军官，结果不幸在一次战斗中被敌军杀死。在他死后的 1924 年，他的算法才得到发表。所以说啊，如果没有第一次世界大战，也许会有更多算法或者定理以楚列斯基名字命名呢。

画外音打断施法：虽然但是，你怎么能出现成的算法？？？

答： 我来解释一下 🤣

1. 首先，这题是可以用大家已学的 OI 知识解决的，就是那个高斯消元的做法。我们也特意把数据改成随机生成的，来让高斯消元可以安全通过 😆 最后确实赛场上有很多人过，所以这说明难度上不成问题。如果你似乎写出了与楚列斯基分解非常接近的算法，那说明你是楚列斯基再世！
2. 其次，我承认确实楚列斯基分解算法现在已经是各大高校的计算机专业在教授《数值分析》这门课会学的算法，大学生来参赛会有一点优势，我们出这题的时候确实也有点担心这一点。但我们并没有出原算法！即使加上前面转化到向量内积的那步，最原始的楚列斯基分解也只能处理点坐标能映射到 $n-1$ 维且不能映射到 $n-2$ 维的情况（即 $A$ 矩阵满秩）。所以下面我们要介绍的是楚列斯基分解算法的一个变种，而非普通的《数值分析》教科书上会介绍的原版算法。
3. 最后，OI 里面特别流行用模意义下的算术来代替实数运算，因为大家往往觉得精度误差问题是“玄学”、“不应该是考察的重点”。然而，升入大学之后，我发现这个理解跟实际有一点点脱节，因为这个时代实在有太多算法是数值算法了。就算是最经典的网络流问题，目前最快的算法也是基于连续优化算法，而非在 Dinic、isap 上继续死磕。所以我觉得，出这道题是个特别好的机会，科普一下在有些情况下，误差是可以分析的，并且是可以通过改进算法来减小误差的（但我意思不是说要在赛场上干这个233）。如果你学有余力，想了解更多，那么你可以找一本《数值分析》教科书看看，不过你需要了解一点点点点线性代数和微积分，所以当然也可以等到大学再学。

下面我们来介绍一下楚列斯基分解算法是如何找到一组满足 (F) 的向量的。与前面的高斯消元不同，这里我们要介绍的楚列斯基分解（的一个变种）不是一个向量一个向量地考虑的，而是一维一维考虑的。

1. 简单情形

首先我们先不考虑询问次数的限制，直接在程序一开始就把整个矩阵 $A$ 给查询出来。同时，我们也不考虑任何数值方面的问题。我们待会儿再来讨论考虑这些因素的话应该怎么做。

第一步，我们考虑第一维。诶，反正什么都没确定过，那就一号向量的方向当做第一维的方向吧！这是不失一般性的 —— 因为如果不是第一维的方向，我们总是可以通过旋转坐标系的方式，把第一维转到该向量的方向上去。

注意到 $A_{1,1}$ 就是一号向量模长的平方，所以一号向量就是 $(\sqrt{A_{1,1}}, 0, \dots, 0)$。这样我们也可以顺便解出 $q = 2, \dots, n-1$ 的每个向量在第一维上的坐标： $$ x_{1,1} = \sqrt{A_{1, 1}} \qquad x_{q,1} = \frac{A_{q,1}}{x_{1,1}} $$ 如此一来，我们似乎就可以删去第一个向量和第一维了。我们构造一个新矩阵 $A'$，满足 $A'_{p,q} = A_{p,q} - x_{p,1} \cdot x_{q,1}$（$2 \le p, q \le n - 1$），然后原问题就从 (F) 转化为了如下问题 (F')： $$ \forall p, q \in \{2, \dots, n -1\}: \qquad \sum_{j=1}^{k} x_{p,j} \cdot x_{q,j} = A'_{p,q}. \tag{F'} $$ 于是递归构造这 $n-2$ 个向量的坐标即可。令递归第 $i$ 层的矩阵是 $A^{(i)}$，那么第 $i$ 层的赋值操作就是： $$ x_{i,i} = \sqrt{A^{(i)}_{i, i}} \qquad x_{q,i} = \frac{A^{(i)}_{q,i}}{x_{i,i}} $$ 递归时构造的新矩阵 $A^{(i+1)}$ 就是 $A^{(i+1)}_{p,q} = A^{(i)}_{p,q} - x_{p,i} \cdot x_{q,i}$（$i+1 \le p, q \le n - 1$）。

每次递归总会少一个向量，所以递归层数不超过 $n$。易见时间复杂度是 $O(n^3)$。

2. 选主元

上面算法有问题吗？有。比如我们没讨论递归的终止条件。除此之外呢？

凡是有除法的地方，我们都需要警惕！注意 $x_{q,i} = A^{(i)}_{q,i}\,/x_{i,i}$ 这里是有个大除法的，那么 $x_{1,1}$ 可不可能是 $0$ 呢？

可能……所以我们得处理这个情况。

处理的方式也很简单，就是仿照高斯消元法里的“选主元”操作，选个模长最大的向量作为主元就好了。

具体来说，就是先扫描一遍所有向量的模长平方 $A^{(i)}_{p,p}$，找到模长最大的那个向量，将其重新编号为 $i$，然后再进行简单情形所描述的算法。

如果最大模长是 $0$（或接近 $0$），那么直接把剩下的维度全部赋值为 $0$ 就好。

易见时间复杂度还是 $O(n^3)$。

3. 减少询问次数

下面我们考虑询问次数有限制的情况。这时我们只要假装往下递归时需要的 $A'$ 已经被我们算出来了就好了，等到最后要执行 $x_{q,i} = A^{(i)}_{q,i}\,/ x_{i,i}$ 时再真正求出 $A^{(i)}$ 的值。

具体来说，我们只用把上面的递归算法变成一个迭代算法。该迭代算法分为若干轮，其中第 $i$ 轮对应于上面递归算法的第 $i$ 层（即用来确定第 $i$ 维的那一层）。 若第 $i$ 层递归时选到的主元是第 $p$ 个向量，那么我们将每次递归时从 $A^{(i)}$ 变到 $A^{(i+1)}$ 的式子逐层展开，就可以知道第 $i$ 层实际执行的操作是： $$ x_{p,i} = \sqrt{A_{p, p} - \sum_{s=1}^{i-1} x_{p,s}^2} \qquad x_{q,i} = \frac{1}{x_{p,i}} \left( A_{q,p} - \sum_{s=1}^{i-1} x_{p,s} \cdot x_{q, s} \right) $$ 其中这里的 $A$ 就是最原始的由 $d$ 定义的 $A$，$q$ 是任意一个还未被选作主元的向量的下标。所以我们只要让这个迭代算法在第 $i$ 轮执行上面的赋值就好了。

选主元的过程是类似的。我们每次会在所有还未被选作主元的向量中，找到 $A_{p, p} - \sum_{s=1}^{i-1} x_{p,s}^2$ 最大的那一个（即上面那个被开根号的东西，也即递归版算法里当前第 $i$ 层的 $A^{(i)}_{p,p}$）

结束条件也是一样的。如果选主元的时候发现 $A_{p, p} - \sum_{s=1}^{i-1} x_{p,s}^2$ 是 $0$（或接近 $0$），就退出。

代码很短，很简洁：https://uoj.ac/submission/543639 （代码里面维度的下标是从 $0$ 到 $k-1$，所以跟上面的式子中的下标有些微小的差异）

当然，你也可以不从递归算法的角度理解。但总之，本质上就是从 $1$ 到 $k$ 按顺序考虑每一维，每次先看看是不是所有的向量都被前面已经构造好的维度完美地表达了。如果是，就退出；如果不是，那么就找一个表达得最不完美的向量，把尚未被表达的部分设为新的维度的方向，然后凑一下看看这个新维度上每个向量的坐标应该是多少，循环往复。

五、误差分析

下面就是激动人心的误差分析环节。

1. 基础知识

回忆一下我们是如何分析时间复杂度的 —— 我们会假设每个操作大概会花费常数的时间，然后数一数操作数，当作时间复杂度。

注意这是个理想模型！实际上不同的运算在计算机上花费的时间是不同的，现代 CPU 也不一定会按顺序执行你的代码。但你算时间复杂度的时候有想这么多吗？木有，因为有些细节完全是可以忽略的，只要粗略地算一个时间复杂度出来就好。

分析精度也是类似。你可能会因为 NOIP 初赛了解一点点浮点型变量的存储规则和运算方法，看起来有点复杂 —— 但你要算精度，就大概带着个大 $O$ 符号算出一个渐进的结果就好。

因为浮点数使用的是科学计数法，所以一切误差都是相对误差。设机器误差是 $\epsilon_0$ （double 大概是 $10^{-15}$），那么任何一个浮点数的四则运算的结果，与准确结果相差的相对误差不会超过 $O(\epsilon_0)$（即绝对误差不超过计算结果的绝对值乘以 $O(\epsilon_0)$）。

举例来说，如果你有两个浮点数 $0 < a, b < R$，那么用浮点数计算 $a+b$ 和直接手算 $a+b$，这两个结果的差不会超过 $O(R \epsilon_0)$；用浮点数计算 $a \cdot b$ 和直接手算 $a \cdot b$，这两个结果的差不会超过 $O(R^2 \epsilon_0)$；

误差是会累积的。如果你把 $n$ 个数不超过 $R$ 的数加起来，误差就会在每一次加法中累积，最后误差就是 $O(nR\epsilon_0)$。

2. 对楚列斯基分解的一个简单分析

好的，那么我们怎么分析楚列斯基分解算法呢？这好像根本不是简单的加减乘除啊……想要计算误差是如何在一步步运算中累积的，似乎需要解一个巨大的递推式，好像算不清楚啊？

这时就要把思路逆转一下：我们不去讨论误差是如何累积的，而是去想，最后算出来的结果满足怎样的性质？

这就好说了。我们再回顾一下 $x$ 是怎么被赋值的： $$ x_{p,i} = \sqrt{A_{p, p} - \sum_{s=1}^{i-1} x_{p,s}^2} \qquad x_{q,i} = \frac{1}{x_{p,i}} \left( A_{q,p} - \sum_{s=1}^{i-1} x_{p,s} \cdot x_{q, s} \right) $$ 但注意啊，上面这个等号是 C++ 的等号，每一个运算都会引入误差，所以数学上是个约等号。

现在我们令 $R$ 为所有 $\lvert x_{p,q} \rvert$ 和所有 $\sqrt{\lvert A_{p,q} \rvert}$ 的最大值。那么就可以这样估计误差了（下面的等号都是真实的数学上的等号）： \begin{align} x_{p,i} &= \sqrt{A_{p, p} - \sum_{s=1}^{i-1} x_{p,s}^2 + O(kR^2 \epsilon_0)} + O(R \epsilon_0) \\ x_{q,i} &= \frac{1}{x_{p,i}} \left( A_{q,p} - \sum_{s=1}^{i-1} x_{p,s} \cdot x_{q, s} + O(k R^2 \epsilon_0) \right) + O(R \epsilon_0) \end{align} 其中最外面的 $O(R \epsilon_0)$ 是因为 <cmath> 的 sqrt 和浮点数除法均保证返回的结果的误差在结果的绝对值乘以 $O(\epsilon_0)$ 以内，而根据 $R$ 的定义，这个结果的绝对值是不超过 $R$ 的。

下面我们看看这组向量的点积和 $A$ 之间到底有多大差距。首先，对 $x_{p,i}$ 平方，有： $$ x_{p,i}^2 = A_{p, p} - \sum_{s=1}^{i-1} x_{p,s}^2 + O(kR^2 \epsilon_0) + 2 \cdot O(R) \cdot O(R \epsilon_0) + O(R^2 \epsilon_0^2) $$ 所以 $\sum_{s=1}^{i} x_{p,s}^2 = A_{p, p} + O(kR^2 \epsilon_0)$。

将 $x_{q,i}$ 乘以 $x_{p,i}$，有： $$ x_{p,i} \cdot x_{q,i} = \left( A_{q,p} - \sum_{s=1}^{i-1} x_{p,s} \cdot x_{q, s} + O(k R^2 \epsilon_0) \right) + O(R^2 \epsilon_0) $$ 所以 $\sum_{s=1}^{i} x_{p,s}^2 = A_{q, p} + O(kR^2 \epsilon_0)$。

最后特殊处理一下迭代结束的情况，我们就可以作出以下结论：对于任意 $1 \le p, q \le n-1$ 都有 $$ \sum_{j=1}^{k} x_{p,j} \cdot x_{q,j} = A_{p,q} + O(k R^2 \epsilon_0), $$ 即 (F) 近似成立，带有 $O(k R^2 \epsilon_0)$ 的绝对误差。

看起来很小，但其实这里有两个问题：

1. 不知道 $R$ 的大小；
2. 精度误差导致 $x_{p,i} = \sqrt{A_{p, p} - \sum_{s=1}^{i-1} x_{p,s}^2}$ 中被开根号的数小于 $0$ 但不约等于 $0$。

当然，如果没有任何误差，$\epsilon_0 = 0$，那么容易看出 $R = O(M)$，其中 $M$ 为所有 $\sqrt{\lvert A_{p,p} \rvert}$ 的最大值，且对负数开根号的情况不会发生。

进一步分析表明，如果 $\epsilon_0$ 足够小，那么 $R = O(M)$ 也成立，且（设置适当的阈值 EPS 的情况下）对负数开根号的情况也不会发生。但把这部分也详细写下来可能篇幅就太长了，感兴趣的同学可以看看这本书的 Theorem 10.14.

所以最后我们就算出了 $O(k M^2 \epsilon_0)$ 的误差上界，似乎很小，所以可以通过本题。

3. 高斯消元会出什么问题？

高斯消元也可以有一个类似的分析，但问题在于，$R$ 的上界可能会很大。

【进入线性代数模式随便扯几句】一般高斯消元的目的是解方程 $Ax = b$，我们关心的是它求出来的近似解 $\hat{x}$ 带回到方程里时的误差是多少，即 $\| A\hat{x} - b \|$ 有多大。要想用类似上面对楚列斯基分解的分析，经过一番分析后会发现，$R$ 至少得是 $A$ 的逆矩阵 $A^{-1}$ 的大小。但是 $A^{-1}$ 可能非常非常大，就会导致最终的精度无法有一个很好的上界。但随机情况的话，$A^{-1}$ 还是挺小的，所以没啥问题。

六、结语

这道题我确实挺喜欢的，既能在 OI 知识范围内解决，也可以引出一些更深的东西。

所以，球面上的距离能表示成某个欧式空间中的距离吗？赶紧拿出你的 AC 程序跑跑看吧 🤣