.

　　在这篇文章中，我们首先来分析一下公式 (4,5) 产生的平面点阵是怎样受转角 $\phi$ 影响的，然后把结论拓展到公式 (1~3) 产生的球面点阵上去。

　　★　 $\begin{align} z_n &= (2n-1)/N - 1 & (1) \\ x_n &= \sqrt{1 - z_n^2} \cdot \cos(2 \pi n \phi) & (2) \\ y_n &= \sqrt{1 - z_n^2} \cdot \sin(2 \pi n \phi) & (3) \end{align}$ ， $\begin{align} x_n &= \sqrt{n} \cos(2\pi n \phi) & (4)\\ y_n &= \sqrt{n} \sin(2\pi n \phi) & (5) \end{align}$

　　分析的过程主要依赖连分式（Continued fraction），在此先列举一些连分式的基础知识。任何一个实数可以展开成如下形式的连分式：

其中各层的分子均为 1； a_0 是的整数部分（取下整）；其它的各个 a_i 都是正整数。如果是有理数，连分式就会像上面一样是有限的；如果是无理数，连分式则是无限的。上面的式子写起来很麻烦，常常简记为 $[a_0; a_1, a_2, \ldots, a_n]$ ，即写成数列的形式， a_i 也称为连分式的项。

　　举几个例子：有理数 10/23 写成连分式是 $0 + \frac{1}{\displaystyle 2 + \frac{1}{\displaystyle 3 + \frac{1}{3}}}$ ，是有限的，简记为 [0; 2, 3, 3] 。而黄金分割比 $\phi = (\sqrt{5} - 1)/2 \approx 0.618$ 写成连分式是 $0 + \frac{1}{\displaystyle 1 + \frac{1}{\displaystyle 1 + \frac{1}{\displaystyle 1 + \frac{1}{\ddots}}}}$ ，是无限的，简记为 $[0; 1, 1, 1, \ldots]$ 。连分式的项不一定有规律，例如圆周率 $\pi$ 的连分式为 $[3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, \ldots]$ （A001203 - OEIS）。

　　把的连分式在某一项处截断，即可得到的一个有理逼近，称为丢番图逼近（Diophantine approximation）。例如，把 $\pi$ 的连分式截断为 [3;7] ，则可以得到「约率」22/7；截断为 [3; 7, 15, 1] ，则可以得到「密率」355/113（请读者自行验证）。一般地，若在连分式中保留到 a_k 这一项，化简后可以得到最简分数 p_k/q_k ，其中分子和分母可以通过以下的递推式得到：

　　★　 $\begin{align} p_k &= a_k p_{k-1} + p_{k-2} & (6) \\ q_k &= a_k q_{k-1} + q_{k-2} & (7) \end{align}$

递推式的初值为 $p_0 = a_0, q_0 = 1; p_{-1} = 1; q_{-1} = 0$ 。

　　丢番图逼近在某种意义下是的最佳逼近。这里「最佳」的意义是：定义 p_k/q_k 的逼近误差 $\epsilon_k = |q_ka-p_k|$ ，则不存在另一个分母不超过 q_k 的分数，逼近误差比 $\epsilon_k$ 小。如果把的各个丢番图逼近排成一个数列 $\frac{p_0}{q_0}, \frac{p_1}{q_1}, \frac{p_2}{q_2}, \ldots$ ，则各项的逼近误差 $\epsilon_0, \epsilon_1, \epsilon_2, \ldots$ 是递减的。

　　（注意，在此我们把误差定义成了 $\epsilon_k = |q_ka-p_k|$ ，而不是 $\left| a - \frac{p_k}{q_k} \right|$ 。前一种定义更适合本文的讨论，因为若把看成每产生一个点旋转的圈数，那么误差 $\epsilon_k$ 的含义就是：连续产生 q_k 个点转过的总圈数，与 p_k 圈差了多少。）

　　有了上面的基础知识，我们来分析一些平面点阵。比如，我们取 $\phi = 0.617$ 。作为准备，我们把 0.617 表示成连分式： [0; 1, 1, 1, 1, 1, 1, 3, 21] ，并计算出它的各阶丢番图逼近及误差：

　　下面是 $\phi = 0.617$ 时，最靠近原点的 400 个点。可以看到，在原点附近，点阵形成螺旋，共有 13 条旋臂；而在外围，点阵呈放射状，共有 47 条射线。不难发现，13 和 47 都是丢番图逼近的分母。我们把螺旋、射线这种有规律的排列称为「模式」（pattern），并用丢番图逼近来命名，例如本图中的两种模式分别称为 8/13 模式和 29/47 模式。

　　200 号点大约处在两种模式的分界线上。观察 200 号点附近的情况：离它最近的点有 153 号、187 号、213 号、247 号，它们的编号与 200 的差恰好也是 13 和 47。把这些点用线连起来，可以发现编号公差为 13 的点列（174 - 187 - 200 - 213 - 226）落在一条 8/13 模式（螺旋）上，而编号公差为 47 的点列（106 -153 - 200 - 247 - 294）落在一条 29/47 模式（射线）上。之所以 8/13 模式产生了明显的旋转，而 29/47 模式不产生明显的旋转，是因为 29/47 的逼近误差远小于 8/13 模式，即：连续产生 47 个点转过的角度与 29 圈的误差，远小于连续产生 13 个点转过的角度与 8 圈的误差。另外还可以发现，在 200 号点附近，8/13 模式上的点的间距随点的编号在增大，而 29/47 模式上的点间距在减小。因此，随着点编号的增大，8/13 模式越来越模糊，而 29/47 模式越来越清晰，导致螺旋被射线所取代。

　　有人会说：「140、166、234、260 这几个点，离 200 号点的距离也不远呀！」这几个点的编号与 200 的差是 34 和 60 —— 发现了没有，它们是 13 与 47 的和与差！把它们对应的 21/34 模式和 37/60 模式也画出来可以发现，这些模式上点的间距比 8/13 模式和 29/47 模式大，因此不如后者显著，可以忽略。这正是因为 21/34 和 37/60 不是 0.617 的最佳逼近 —— 21/34 的逼近误差大于 8/13，37/60 的逼近误差大于 29/47。

　　经过上面的分析，我们可以得到如下的经验：

每个最佳逼近，对应着一种比较显著的模式；
模式上的点间距越小，模式越显著；
同一个模式上的点间距随位置会发生变化，导致最显著模式的切换。

为了发掘模式切换的规律，下面我们来定量地计算模式上点的间距。

　　考虑号点，以及它所处的 p_k/q_k 模式。在该模式上，下一个点的编号是 n+q_k 。把这两个点的编号代入公式 (4,5)，就能得到两点的坐标和距离了。但是这样得到的公式太复杂，我们换一个角度来思考。号点距原点的距离为 $r=\sqrt{n}$ ，在号点附近，编号每增加 1，距原点的距离增加 $\frac{\mathrm{d}r}{\mathrm{d}n} = \frac{1}{2\sqrt{n}} = \frac{1}{2r}$ 。当 $q_k \ll n$ 时，从号点到 n+q_k 号点，径向距离就是 $\frac{q_k}{2r}$ 。而从号点到 n+q_k 号点转过的总圈数为 $q_k\phi$ ，它大约是 p_k 整圈，误差为 $\epsilon_k = |q_k \phi - p_k|$ ，于是这两点间的切向距离就是 $2\pi r \epsilon_k$ 。当 $\epsilon_k \ll 1$ 时，可以把这个切向距离看成与径向垂直的直线段，由此得到 p_k/q_k 模式在距原点处的点间距：

　　★　 $\begin{align} d_k &= \sqrt{\left(\frac{q_k}{2r}\right)^2 + (2\pi r \epsilon_k)^2} & (8) \end{align}$

这里我选择距原点距离而不是点的编号为自变量，显得更直观。

　　如果选定某一种模式（即固定 $q_k, \epsilon_k$ ），则 d_k 是的函数。根号下的第一项随递减，第二项随递增，整体上则呈现先减后增的「V」字形。 d_k 的最小值为 $d_{k,\min} = \sqrt{2 \cdot \frac{q_k}{2r} \cdot 2\pi r \epsilon_k} = \sqrt{2\pi q_k \epsilon_k}$ ，它在 $\frac{q_k}{2r} = 2\pi r \epsilon_k$ 即 $r_{k,\min} = \sqrt{\frac{q_k}{4\pi\epsilon_k}}$ 处取得。由此可知，从原点向外，各种模式会先变得越来越显著，这是因为半径的增长越来越慢，使得模式上的点靠得越来越近；然后会变得越来越模糊，因为转角 $q_k\phi$ 与 p_k 整圈的微小偏差被越来越大的半径放大成越来越远的距离。这个规律也可以通俗地描述成：每个模式先是形成越来越清晰的射线，然后射线弯曲成螺旋，最后「被风吹散」。根据 $r_{k,\min} = \sqrt{\frac{q_k}{4\pi\epsilon_k}}$ ，值越大的模式，分母 q_k 越大，而误差 $\epsilon_k$ 越小，故使得该模式最显著的半径 $r_{k,\min}$ 越大。也就是说，值越大的模式，就会在离原点越远处显现，这造成了模式的交替。

　　我们用上面 $\phi = 0.617$ 时的平面点阵检验一下上面的结论。下表计算了各个模式最显著处的半径 $r_{k,\min}$ （称为「最显著半径」）和此处的点间距 $d_{k,\min}$ （称为「最密点间距」）。 $k \le 5$ 的各个模式（包括 0/1, 1/1, 1/2, 2/3, 3/5, 5/8），误差都太大，同时 $r_{k,\min}$ 太小、 $d_{k,\min}$ 太大，所以在图上显现不出来。8/13 模式在半径约为 7 处（50 号点附近）最显著，这与图示相符；29/47 模式在半径约为 61 处最显著，但前文的图中只画出了半径为 20 以内的部分，所以这一点就看不出来了。从表中可以看到 8/13 模式的最密点间距为 1.310，而 29/47 模式的最密点间距为 0.543，所以若把点阵继续画下去， 29/47 模式会变得比 8/13 模式更显著。最后还有一个模式 617/1000，它将在 29/47 模式消散后显现出来。并且因为 617/1000 就是 $\phi$ 的精确值，此模式将永不消散，并越来越显著。

　　令 $d_k = d_{k+1}$ ，可以求出模式 p_k/q_k 向 $p_{k+1}/q_{k+1}$ 过渡处的半径。仍然研究 $\phi = 0.617$ 的例子，取 k=6 ，可以求得模式 8/13 与 29/47 之间的过渡发生在半径为 13.090 处（170 号点附近）。我们前面观察过的 200 号点离这儿也不远，从图上可以看到这里确实发生了模式的过渡。同样，取 k=7 ，可以求得模式 29/47 到 617/1000 的过渡发生在半径为 281.939 处，也就是说，要画到 80000 号点，才能观察到这个过渡。

　　为了让大家一饱眼福，我把点阵画到了 30 万个点，并截取了横轴正半轴周围的局部。请点开大图观察，否则图片上会产生摩尔纹（Moiré pattern）。从图中可以清楚地看到原点附近的 8/13 模式（螺旋），它在半径为 13 处过渡为 29/47 模式，后者在半径 60 左右最显著，且比 8/13 模式更显著。29/47 模式的竖纹在半径 250 处基本消散，到半径 350 处开始产生 617/1000 模式的横纹。

　　当 $\phi = 0.617$ 时各模式的点间距随半径的变化可以总结成下图，注意横、纵轴均为对数刻度。从图中可以清楚地看出 8/13、29/47、617/1000 三种主要模式的最显著半径、最密点间距，以及模式之间过渡的位置。一图胜千言，下文中我将主要使用这种图（称为模式图）来说明模式的变化。

　　下面，我们来探究连分式中项的大小对模式的影响。我们选取 $\phi = \pi - 3 \approx 0.1416$ 。这个值的连分式为 $[0; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, \ldots]$ ，其中 15、292 这样的大项，会带给我们新的发现。

　　先画个点阵看看：

图中可以观察到两种模式：中心的模式有 7 条腿，外围的模式腿太多了，数不清。

　　按上文给出的公式，计算 $\phi$ 的丢番图逼近序列、各个逼近的误差，以及各个模式的最显著半径和最密点间距，并作出模式图，如下所示。

　　从模式图中观察到，1/7 和 16/113 这两个模式十分显著，它们在半径为 30 处左右发生过渡，这与点阵图完全相符。注意，在这两个模式发生过渡的时候，有一个 15/106 模式被完全掩盖过去了。16/113 这个模式会一直持续到半径接近 10000 的时候才会过渡到 4703/33215 这个模式，哪怕画大图也无缘查看了。

　　1/7 和 16/113 这两个模式，正是连分式中 15 和 292 这两个大项产生的。这里的「大」，其实并不是指绝对数值很大，而是指相对于周围的项来说很大。回忆模式最密点间距公式 $d_{k,\min} = \sqrt{2\pi q_k \epsilon_k}$ ，其中 q_k 是递增的， $\epsilon_k$ 是递减的。再回忆丢番图逼近中分母的递推式 $q_k = a_k q_{k-1} + q_{k-2}$ ：当 a_k 是一个大项时， q_k 与 $q_{k-1}$ 相比，就会猛然增大。但当 a_k 很大时，若把连分式截断至 $a_{k-1}$ ，那么舍弃的 $\frac{1}{a_k + \ddots}$ 就很小，所以 $\epsilon_{k-1}$ 与之前的 $\epsilon_{k-2}$ 相比就会猛然减小。误差的减小比分母的增大早一步发生，这就使得 $d_{k-1, \min}$ 明显小于相邻模式的最密点间距，产生一个十分显著的模式 $p_{k-1}/q_{k-1}$ 。另外观察模式图可以发现，各个模式的「V」字形的两臂斜率都是一样的，所以越是显著的模式，覆盖的半径范围往往也越大。

　　上面的计算告诉我们，连分式中每出现一个大项，在它之前截断就会产生一个十分显著的模式。那么，若要不产生显著的模式，连分式中就不能出现大项，所以令所有的 a_k 都等于 1 就是最好的选择了。这就给出了最优的转角 —— 黄金角。当 $\phi = (\sqrt{5} - 1) / 2 \approx 0.618$ 时，模式图如下（注意，各模式的分子、分母都是菲波那契数）：

　　可以看出，各个模式的显著程度都差不多，并且模式切换频繁。观察下面点阵图的中心，似乎能看出许多层层叠叠的花瓣，这正是模式频繁切换的体现。到了外围，模式持续的长度就长了些（注意模式图的横坐标是对数刻度），能看出一些模式了；这些模式都呈螺旋而不是射线状，因为各阶丢番图逼近的误差都不小。

　　模式切换频繁正是点阵混乱和密集的原因：混乱是因为任何一种模式（规律）都不会持续很久；密集是因为在任一半径处，都并没有单独一种模式的点间距比别的模式小得多，所以不会出现大的缝隙。

　　可以证明，对于黄金分割比 $\phi$ ，丢番图逼近 p_k/q_k 的误差 $\epsilon_k \approx \frac{1}{\sqrt{5}q_k}$ ，于是各模式的最密点间距就都约等于 $\sqrt{2 \pi / \sqrt{5}} \approx 1.676$ 。这其实是用形如 (4,5) 的公式能生成的点阵中，模式最密点间距的「上界」—— 加引号是因为，可能出现个别模式的最密点间距大于这个值，但这样的个别模式至多只有有限个。

　　话说回来，公式 (4,5) 并不是平面上点阵密铺问题的最优解。在保证同样点阵密度（每 $\pi$ 面积上有一个点）的情况下，正方形网格的最密点间距是 $d_\text{square} = \sqrt{\pi} \approx 1.772$ ，而最优解六边形网格的最密点间距是 $d_\text{hex} = \sqrt{2\pi/\sqrt{3}} \approx 1.905$ （如下图）。这两种网格满足均匀、密集，但不满足混乱，因为同一种排列规律遍布了整个平面。从 1.905 到 1.676，就是「混乱」需要付出的代价。

　　现在来看一下我在试探过程中发现的另一个不错的 $\phi$ 值： $\phi = \sqrt{2} - 1 \approx 0.414$ 。此时的点阵图与模式图如下：

　　把 $\phi = \sqrt{2} - 1 \approx 0.414$ 写成连分式，是 $[0; 2, 2, 2, \ldots]$ ，各项都相等，所以没有「大项」。模式图上也显示出，各个模式的最密点间距都相等，且模式切换频繁。但是，与黄金分割比相比， $\phi = \sqrt{2} - 1 \approx 0.414$ 的各个模式的最密点间距偏小（约为 $\sqrt{\pi/\sqrt{2}} \approx 1.490$ ），故模式还是稍微明显了一些；模式之间的切换也不如黄金分割比频繁（注意两张模式图的横、纵坐标范围是相同的）。黄金分割比以微弱的优势战胜 $\sqrt{2} - 1$ 。

　　上面我们详细地分析了 $\phi = 0.617,\,\, \pi - 3,\,\, \frac{\sqrt{5} - 1}{2},\,\, \sqrt{2} - 1$ 各值时的点阵图与模式图，得到了如下一些结论：

每个点转过的圈数 $\phi$ 的每个丢番图逼近对应着一种模式；
随半径增加，每种模式先变得显著，然后消散，存在一个最显著半径；
分母越大的模式，最显著半径也越大，所以各模式会依次显现；
当连分式中有大项的时候，大项之前的那个丢番图逼近对应的模式会非常显著，且持续的长度也很长；
黄金分割比的连分式中所有项都是 1，这导致没有显著的模式且模式切换频繁，于是点阵显得密集而混乱。

　　终于到了把结论拓展到球面的时候了。在平面上，我们选择了到原点的距离作为自变量；在球面上，我们选择纬度 $\theta$ 作为自变量。设整个点阵的点数为，考虑模式 p_k/q_k 上的号点和 n+q_k 号点。在纬度为 $\theta$ 处，这两个点在经线方向上的距离为 $\frac{2q_k}{N\cos\theta}$ ，在纬线方向上的距离为 $2 \pi \epsilon_k \cos\theta$ 。由此可得 p_k/q_k 模式在竖坐标为处的点间距：

　　★　 $\begin{align} d_k &= \sqrt{\left(\frac{2q_k}{N \cos\theta}\right)^2 + (2\pi \epsilon_k \cos\theta)^2} & (9) \end{align}$

此式的最小值为 $d_{k,\min} = \sqrt{8\pi q_k \epsilon_k / N}$ ，在 $\theta_{k,\min} = \pm \arccos\sqrt{\frac{q_k}{N\pi\epsilon_k}}$ 处取得。作为对比，平面上 p_k/q_k 模式的最小点间距为 $d_{k,\min} = \sqrt{2\pi q_k \epsilon_k}$ ，在 $r_{k,\min} = \sqrt{\frac{q_k}{4\pi\epsilon_k}}$ 处取得。两个最小点间距的公式是相似的，只是系数不同；但最显著位置的公式则有形式上的差别：球面上的公式多了个反余弦。反余弦函数是单调递减的，且定义域有限，这说明，各个模式是从两极到赤道依次显现的，而有些高阶的模式，哪怕到了赤道处也来不及显现。除此以外，平面点阵的各个结论都适用于球面点阵。

　　例如，下面是 $\phi = 0.617$ 时，在球面上取 1000 个点生成的点阵，以及相应的模式图。可以看到 8/13 模式（螺旋）在纬度 64 度时最显著，并在 35 度左右过渡成 29/47 模式（射线）。但在低纬度区域，8/13 模式并没有消散得很厉害，它与 29/47 模式共同组成了方形网格。617/1000 模式来不及显现出来。

　　下图是 $\phi$ 取黄金分割比时，由 1000 个点组成的球面点阵和模式图。与平面的情况类似，这里没有特别显著的模式，且模式切换频繁，所以点阵的分布密集且混乱。不过，由于 arccos 定义域的影响，21/34 这个模式从北纬 40 度一直到南纬 40 度都是最显著的；但与此同时，相邻的 13/21 模式和 34/55 模式的显著程度也差不多，所以在低纬度区域的点阵中能找到正方形或六边形网格的影子。

ITPub博客

10561 怎样在球面上「均匀」排列许多点？（下）

最新文章