朝花夕拾

轻飘飘的旧时光就这么溜走 转头回去看看时已匆匆数年

高颜值小程序:用 Python + POV-Ray 绘制多面体和多胞体

2018-05-21


本文要介绍的是我写的另一个高颜值的、逼格满满的、脱离了低级趣味的小程序:用 Python 和 POV-Ray 绘制各种三维多面体和四维多胞体,代码仍然放在 github 上的 pywonderland 项目中。

这个项目目前还是一个半成品,但是已经可以做很多有意思的事情,下面是用这个程序渲染的一些例子:

例子

  1. 5 种 Platonic 正多面体

  2. 13 种 Archimedean 均匀多面体

  3. 棱柱和反棱柱:

  4. 6 种正多胞体

  5. 将以上 6 种正多胞体通过截断得到的均匀多胞体:(由于截断的方式较多,仅对每个正多胞体列出其一种截断)

  6. 四维棱柱和双棱柱:

  7. 4 种 Kepler-Poinsot 正多面体:

我还制作了一个视频,尝试重现精彩的数学视频 Dimensions 中的效果:下面视频展示的是一个 4 维的多胞体 120-cell,它有 600 个顶点,1200 条边,720 个面,120 个胞腔。所有这些顶点、边、面、胞腔在 4 维空间中都是一样大小的,只不过在球极投影到 3 维空间时发生了形变。

这些图都是在 Python 中做好计算以后,将多面体/多胞体的顶点、边、面的数据导出到 POV-Ray 中渲染得到的。你也完全可以通过自己改写代码中的 POV-Ray 的部分来渲染得出不同的效果,当然前提是需要掌握 POV-Ray 的场景描述语言,不过这属于另一段故事,不在本文的讨论范围内。

下面介绍程序背后的数学原理。

这些图画的都是什么?

严格的讲,这些图都是三维或者四维欧式空间凸/非凸均匀多面体/多胞体。其中三维的情形叫做多面体,四维的情形叫做多胞体。这里有几个关键词需要注意:凸/非凸、欧式空间、均匀。

"凸" 比较好理解,就是指多面体/多胞体上任意两点间的连线仍然属于此多面体/多胞体或者其内部。有 "凸" 当然就有 "非凸",典型的如 small stellated dodecahedron:(下面这个图是用本程序渲染的!)

在三维空间中一共有 18 种凸的均匀多面体 (不含两个无穷大类棱柱和反棱柱),57 种非凸的均匀多面体。这个程序目前只能画凸多面体和部分非凸多面体,看起来要使得它能够画所有的 3D/4D 均匀多面体还有不少活儿要做,暂时放到以后吧。

"欧式空间" 也很好理解,这里特别强调它是因为在三维和四维实向量空间中还可以定义其它度量 (比如双曲度量),这时空间中的度量是 "弯曲" 的,从而其中的 "正多面体" 会呈现另外一种样子,最著名的例子是数学软件 mathemetica 的 logo "Spikey",它是一个三维双曲空间中的正十二面体:

"均匀" 这个词就不太好理解了。简单说,就是所有顶点都一样,且每个面都是正多边形 (对多胞体还要加上每个胞腔都是均匀多面体)。要准确解释什么叫 "所有顶点都一样",就要用到群论的语言:一个多面体/多胞体 \(P\) 的对称群 \(G\) 是欧式空间中一组正交变换构成的有限群,且 \(G\) 作用在 \(P\) 上保持 \(P\) 不变。"所有顶点都一样" 的严格说法是 "\(G\) 传递地作用在 \(P\) 的顶点集上",即对 \(P\) 的任何两个顶点 \(u,v\),都存在 \(g\in G\)\(g\) 把顶点 \(u\) 映射为顶点 \(v\)

这些图是怎么画出来的?

前面说过,这些多面体/多胞体的数据是在 Python 中计算好以后,导出到 POV-Ray 中渲染得到的。这其中的重点在 Python 的部分,POV-Ray 部分的计算虽然也略复杂,但是用到的知识不超出高中几何的范围,没有很大的难度。

这些多面体/多胞体看起来样子大不相同,但实际上它们都可以用同一种方法计算出来,叫做 Wythonff 构造法,也称 "万花筒构造法"。这个构造法的原理跟我们小时候玩的万花筒的原理是一样的:即在空间中放置若干过原点的反射平面 (镜子),这些镜面之间的夹角是精心设计好的,都形如 \(\pi-\pi/p\),其中 \(p\) 为有理数。在空间中选定一个初始顶点 \(v_0\),将 \(v_0\) 关于这些镜子反复作反射变换,得到的全部镜像就是多面体/多胞体的顶点。如果某个镜像点 \(v\) 关于第 \(i\) 面镜子反射后得到的镜像是 \(v'\),则 \((v,v')\) 构成一条边,我们把它染成 \(i\) 号色。如果一个镜像点 \(v\) 先关于镜面 \(i\) 作反射,再关于镜面 \(j\) 作反射,则由于两个反射变换的复合是一个旋转变换,\(v\) 实际上是绕着某个面的中心和原点的连线作了一次旋转,旋转的角度为 \(2\pi/m\) (假设镜面 \(i\) 和镜面 \(j\) 的法向量夹角是 \(\pi-\pi/m\)),重复此旋转 \(m\) 次即可得到多面体/多胞体的一个面。

这里的关键之处有两个:

  1. 对于不同的凸均匀多面体/多胞体,应该如何放置这些镜面和选择初始顶点?
  2. 摆好镜面和初始顶点以后,怎样求出初始顶点的所有镜像?

第一个问题的答案叫做 Coxeter-Dynkin 图,Coxeter-Dynkin 图是一个带标记信息的无向图,它编码了多面体/多胞体的全部信息。这里 "编码了全部信息" 的意思是只要知道了多面体/多胞体对应的 Coxeter-Dynkin 图,就可以求出该多面体/多胞体的所有数据 (仅缩放大小和在空间中的摆放位置除外)。每个均匀多面体/多胞体都有一个 Coxeter-Dynkin 图与之对应,但是不同的 Coxeter-Dynkin 图可能描述的是相同的多面体/多胞体。

比如正方体的 Coxeter-Dynkin 图为:

我们来解释这个图的含义:

在一个 Coxeter-Dynkin 图中,每个顶点代表一面镜子,在上图中有三个顶点,所以有三面镜子。将这三面镜子从左到右依次记作 \(m_0, m_1, m_2\),顶点之间的边记录了镜子间的夹角:

  1. 若两个镜面之间的夹角为 \(\pi/2\) 则它们之间没有边相连。
  2. 若两个镜面之间的夹角为 \(\pi-\pi/3\) 则它们之间用一条无标记的边相连。
  3. 若两个镜面之间的夹角为 \(\pi-\pi/m, m>3\) 则它们之间用一条标号为 \(m\) 的边相连。

此外用 "圈出" 的镜面来标记初始顶点的位置,一个镜面被圈出当且仅当初始顶点不在这个镜面上

从而在正方形的情形 \(\langle m_0,m_1\rangle=\pi-\pi/4\)\(\langle m_1,m_2\rangle=\pi-\pi/3\)\(\langle m_0,m_2\rangle=\pi/2\)。初始顶点落在 \(m_1\)\(m_2\) 上但是不属于 \(m_0\)

于是这三面镜子可以这样摆放:

  1. 镜子 \(m_0\) 的法向量可以随意,比如 \(n_0=(1, 0, 0)\)
  2. 镜子 \(m_1\) 的法向量 \(n_1\)\(n_0\) 夹角为 \(3\pi/4\),于是 \(n_1\) 可以取为 \((\cos\dfrac{3\pi}{4}, \sin\dfrac{3\pi}{4}, 0)\)
  3. 镜子 \(m_2\) 的法向量 \(n_2\)\(n_0\) 垂直,所以 \(n_2\) 形如 \((0,y_3,z_3)\),它与 \(n_1\) 夹角是 \(2\pi/3\),所以 \(y_3 \sin\dfrac{3\pi}{4}=\cos\dfrac{2\pi}{3}\)\(z_3=\sqrt{1-y_3^2}\),解出 \(y_3, z_3\) 即可。

要选择一个落在 \(m_1\)\(m_2\) 上但是不落在 \(m_0\) 上的初始点 \(v_0\),我们可以让 \(v_0\) 到平面 \(m_1\)\(m_2\) 的距离为 0,到平面 \(m_0\) 的距离为 1,即 \[\begin{align*}\langle v_0, n_0\rangle=1,\\\langle v_0, n_1\rangle=0,\\\langle v_0, n_2\rangle=0.\\\end{align*}\]

求解这个线性方程组即可。

再比如 omnitruncated dodecahedron 的 Coxeter-Dynkin 图为

即初始顶点不属于三个镜面中的任何一个。

我们前面提到过,要使得初始顶点的所有镜像恰好构成一个均匀多面体/多胞体,镜子之间的夹角必须精心设置,这实际上只有有限种可能。换句话说,只有有限几种 Coxeter-Dynkin 图可以给出均匀多面体/多胞体,它们分别对应 rank 为 3 或者 4 的几种有限 Coxeter 群。在维基百科上完整的列出了每种均匀多面体/多胞体对应的 Coxeter-Dynkin 图,这里就不再专门列举了,但是特别指出两点:

  1. Coxeter-Dynkin 图的标号完全决定了多面体的对称性,而初始顶点的位置则决定了多面体的截断类型。
  2. 对偶的多面体具有相同的 Coxeter-Dynkin 图,只不过要把边的标号从右到左反过来。比如正八面体和正方体的 Coxeter-Dynkin 图是一样的,但是边的标号是 (3, 4)。

第二个问题的答案叫做 Todd-Coxeter 算法,展开讲的话比较复杂,我们单列一节专门谈谈。

有限表现群和 Todd-Coxeter 算法

怎样求出初始顶点在所有镜子中的镜像?答案看起来好像很显然:只要反复地将初始顶点关于每个镜面作反射,算出得到的镜像点的坐标,并与之前得到的点的坐标相比较 (浮点数向量比较相等需要事先规定在一定的误差范围内考虑),直到不再有新的镜像点出现为止,得到的全部镜像点不就构成顶点集了吗?这个方法确实可行,但是既笨又丑陋:它完全没有用到多面体/多胞体具有对称性这一事实!

不过客观讲,这种无脑反射的方法也有其价值,比如用 glsl 来渲染时这是一个很自然的选择,这个 shadertoy 上的项目就是采用的这种途径,效果也很棒。

这个程序采用的是另外一种纯群论的 "符号计算" 的途径,这个方法可以 "精准" 地得出任何一个顶点/边/面的所有信息,代价就是涉及的数学略复杂,需要读者学过抽象代数中群论的知识。在介绍这个途径之前我们回忆一下关于群在集合上的作用的轨道 - 稳定化子定理:

轨道 - 稳定化子定理:设群 \(G\) 传递地作用在集合 \(S\) 上,\(x\in S\)\(x\) 的稳定化子群是 \(H\)。则集合 \(S\)\(G/H\) 中的右陪集之间有一一对应:\(Hg\to x\cdot g\)

注意:和一般的约定不同,这里群在集合上的作用为 "作用在右边",主要是为了编程方便,实际上左边右边都一样。

这个定理告诉我们,如果群 \(G\) 传递地作用在一个集合 \(S\) 上,而且对 \(S\) 中某个元素 \(x\) 我们知道了它的稳定化子群 \(H\),则只要对 \(G/H\) 的每个陪集代表元 \(g\),将 \(g\) 作用在 \(x\) 上就可以得到整个 \(S\)

于是给定一个均匀多面体/多胞体 \(P\),要求出其全部顶点集合,我们只要:

  1. \(P\) 的 Coxeter-Dynkin 图出发定出 \(P\) 的对称群 \(G\) 和初始顶点 \(v_0\)
  2. 定出 \(v_0\) 的稳定化子群 \(H\) 并求出 \(G/H\) 的一组陪集代表元;
  3. \(G/H\) 中的每个陪集代表元作用在 \(v_0\) 上即得 \(P\) 的全部顶点。

我们仍然以正方体为例来讲解:回忆正方体的 Coxeter-Dynkin 图为

仍然记三个镜面为 \(m_0,m_1,m_2\),其法向量分别为 \(n_0,n_1,n_2\)\(\rho_0,\rho_1,\rho_2\) 分别是关于它们的反射变换,其中 \(\rho_i\) 对应的矩阵为 \(M_i=I - 2n_in_i^T\) (见 Householder 变换)。

正方体的对称群 \(G\)\(\rho_0,\rho_1,\rho_2\) 这三个基本反射生成,其表现为: \[G = \langle\rho_0,\rho_1,\rho_2\ |\ \rho_0^2=\rho_1^2=\rho_2^2=(\rho_0\rho_1)^4=(\rho_1\rho_2)^3=(\rho_0\rho_2)^2=1\rangle.\] 这是因为由于 \(\rho_0, \rho_1\) 是两个夹角为 \(3\pi/4\) 的反射,所以 \(\rho_0\rho_1\) 是一个角度为 \(3\pi/2\) 的旋转,旋转轴为 \(m_0\)\(m_1\) 的交线,从而 \((\rho_0\rho_1)^4=1\)\(\rho_1\rho_2,\rho_0\rho_2\) 的情形是类似的。

利用 Todd-Coxeter 算法 (后面有解释) 不难求出这个群包含 48 个元素,罗列如下: \[\begin{array}{lll}e&\rho_{0}&\rho_{0}\rho_{1}\\\rho_{0}\rho_{1}\rho_{0}&\rho_{0}\rho_{1}\rho_{0}\rho_{1}&\rho_{1}\rho_{0}\rho_{1}\\\rho_{1}\rho_{0}&\rho_{1}&\rho_{0}\rho_{2}\\\rho_{2}&\rho_{1}\rho_{2}&\rho_{1}\rho_{2}\rho_{1}\\\rho_{2}\rho_{1}&\rho_{0}\rho_{1}\rho_{2}&\rho_{0}\rho_{1}\rho_{2}\rho_{1}\\\rho_{0}\rho_{2}\rho_{1}&\rho_{0}\rho_{1}\rho_{0}\rho_{2}&\rho_{0}\rho_{1}\rho_{0}\rho_{1}\rho_{2}\\\rho_{0}\rho_{1}\rho_{0}\rho_{1}\rho_{2}\rho_{1}&\rho_{0}\rho_{1}\rho_{0}\rho_{2}\rho_{1}&\rho_{1}\rho_{0}\rho_{1}\rho_{2}\\\rho_{1}\rho_{0}\rho_{2}&\rho_{1}\rho_{0}\rho_{2}\rho_{1}&\rho_{1}\rho_{0}\rho_{1}\rho_{2}\rho_{1}\\\rho_{2}\rho_{1}\rho_{0}&\rho_{2}\rho_{1}\rho_{0}\rho_{1}&\rho_{0}\rho_{2}\rho_{1}\rho_{0}\rho_{1}\\\rho_{0}\rho_{2}\rho_{1}\rho_{0}&\rho_{1}\rho_{0}\rho_{2}\rho_{1}\rho_{0}&\rho_{1}\rho_{0}\rho_{2}\rho_{1}\rho_{0}\rho_{1}\\\rho_{1}\rho_{2}\rho_{1}\rho_{0}\rho_{1}&\rho_{1}\rho_{2}\rho_{1}\rho_{0}&\rho_{0}\rho_{1}\rho_{0}\rho_{2}\rho_{1}\rho_{0}\\\rho_{0}\rho_{1}\rho_{0}\rho_{2}\rho_{1}\rho_{0}\rho_{1}&\rho_{0}\rho_{1}\rho_{2}\rho_{1}\rho_{0}\rho_{1}&\rho_{0}\rho_{1}\rho_{2}\rho_{1}\rho_{0}\\\rho_{1}\rho_{0}\rho_{1}\rho_{2}\rho_{1}\rho_{0}&\rho_{1}\rho_{0}\rho_{1}\rho_{2}\rho_{1}\rho_{0}\rho_{1}&\rho_{0}\rho_{1}\rho_{0}\rho_{1}\rho_{2}\rho_{1}\rho_{0}\rho_{1}\\\rho_{0}\rho_{1}\rho_{0}\rho_{1}\rho_{2}\rho_{1}\rho_{0}&\rho_{2}\rho_{1}\rho_{0}\rho_{1}\rho_{2}&\rho_{1}\rho_{2}\rho_{1}\rho_{0}\rho_{1}\rho_{2}\\\rho_{0}\rho_{2}\rho_{1}\rho_{0}\rho_{1}\rho_{2}&\rho_{0}\rho_{1}\rho_{2}\rho_{1}\rho_{0}\rho_{1}\rho_{2}&\rho_{1}\rho_{0}\rho_{2}\rho_{1}\rho_{0}\rho_{1}\rho_{2}\\\rho_{1}\rho_{0}\rho_{1}\rho_{2}\rho_{1}\rho_{0}\rho_{1}\rho_{2}&\rho_{0}\rho_{1}\rho_{0}\rho_{2}\rho_{1}\rho_{0}\rho_{1}\rho_{2}&\rho_{0}\rho_{1}\rho_{0}\rho_{1}\rho_{2}\rho_{1}\rho_{0}\rho_{1}\rho_{2}\end{array} \]

由于在正方形的 Coxeter-Dynkin 图中只有镜面 \(m_0\) 是被圈出的,即初始顶点 \(v_0\) 落在 \(m_1\)\(m_2\) 上,但不属于 \(m_0\),所以镜面反射 \(\rho_1,\rho_2\) 作用在 \(v_0\) 上保持其不动,\(\rho_0\)\(v_0\) 映射为其关于 \(m_0\) 的对称点,所以 \[H=\langle \rho_1, \rho_2\ |\ \rho_1^2=\rho_2^2=(\rho_1\rho_2)^3=e\rangle\] 就是 \(v_0\) 的稳定化子群 1。显然 \(H\) 恰好就是二面体群 \(D_3\) (即正三角形的对称群),所以 \(|H|=6\),从而 \(|G/H|=8\)。 仍然利用 Todd-Coxeter 算法可得其一组右陪集代表元为 \[\begin{array}{llll}e&\rho_{0}&\rho_{0}\rho_{1}&\rho_{0}\rho_{1}\rho_{0}\\\rho_{0}\rho_{1}\rho_{2}&\rho_{0}\rho_{1}\rho_{0}\rho_{2}&\rho_{0}\rho_{1}\rho_{0}\rho_{2}\rho_{1}&\rho_{0}\rho_{1}\rho_{0}\rho_{2}\rho_{1}\rho_{0}\end{array} \] 将它们作用在 \(v_0\) 上即可得到正方体的 8 个顶点。例如 \(\rho_0\rho_1\) 作用在 \(v_0\) 上为 \[v_0(\rho_0\rho_1)=(v_0\rho_0)\rho_1=(v_0M_0)\rho_1=v_0M_0M_1.\] 其中 \(v_0\) 写成行向量的形式, 每个 \(M_i\) 是对称矩阵。

注意:这里是把每个生成元从左到右依次作用在 \(v_0\) 上。

求边和面的过程是类似的,比如我们想求出所有关于第 \(i(i=0,1,2)\) 个镜面反射生成的类型为 \(i\) 的边,可以这样做:

  1. 检查初始顶点 \(v_0\) 是否落在镜面 \(i\) 上。如果答案为是,那么关于此镜面的反射保持 \(v_0\) 不变,此多面体不含类型 \(i\) 的边。如果答案是否,设 \(v_0\) 关于此镜面的镜像为 \(v_1\),则 \((v_0, v_1)\) 构成一条类型为 \(i\) 的边 \(e\)
  2. 关于镜面 \(i\) 的反射 \(\rho_i\)\(v_0\)\(v_1\) 互换,从而保持 \(e\) 不变,因此是 \(e\) 的稳定化子群 \(H\) 的生成元:\(H=\langle \rho_i \rangle\) 2
  3. 求出 \(G/H\) 的一组陪集代表元并作用在 \(e\) 上得出全部类型为 \(i\) 的边。(需要剔除形如 \((u, v), (v,u)\) 这种重复的边)

求面的情形复杂一些,基本原理是这样的:对 \(i\ne j\),如果初始顶点 \(v_0\) 不同时属于镜面 \(i\) 和镜面 \(j\),则旋转 \(\rho_i\rho_j\) 就可以生成一个面 \(f\)\(f\) 在这个旋转下不变。当然这里面有一些细节需要考虑,有时 \(f\) 的稳定化子群 \(H=\langle \rho_i\rho_j \rangle\),有时 \(H=\langle \rho_i,\rho_j \rangle\) 3。仍然由此计算 \(G/H\) 的一组陪集代表元,将每个代表元作用在 \(f\) 上即得所有对应旋转 \(\rho_i\rho_j\) 的面。

总之关键的步骤都是给定群 \(G\) 和某个子群 \(H\),求 \(G/H\) 的一组陪集代表元

这里用到的算法叫做 Todd-Coxeter 算法

Todd-Coxeter 算法在许多抽象代数或者群论的教材都有,用到的数学知识并不复杂,但是完整描述并证明一份程序实现的细节还是很费功夫的,恐怕要好几页纸才能说清楚。限于篇幅,我这里仅用正方体的情形为例说明算法的步骤,具体的证明和更多的细节读者请参考

Handbook of Computational Group Theory, Holt, D., Eick, B., O'Brien, E.

中的 coset enumeration 一章。这是目前讲解 Todd-Coxeter 算法最详尽的文献。

Todd-Coxete 算法非常类似玩数独游戏,不过这里我们要填的表是一个变化的二维数组 \(T\)\(T\) 的行 \(i\) 代表第 \(i\) 个右陪集,\(T\) 的列 \(j\) 代表第 \(j\) 个生成元 \(\rho_j\)\(T[i][j]\) 的值等于 \(\rho_j\) 右乘以第 \(i\) 个陪集后得到的陪集。初始时,我们知道肯定有一个陪集,就是 \(H\) 自身,还有没有其它的陪集我们不清楚。算法的 "终极目的" 就是根据 \(G\) 的生成关系和 \(H\) 的生成字这些信息来发现新的陪集并填入表中,直到无法找到新的陪集为止。这时得到的 \(T\) 实际上是 \(G/H\) 的 schreier 图的邻接矩阵,它记录了 \(G/H\) 中的陪集之间的乘法关系,由 \(T\) 出发我们很容易求出这些陪集的 word 表示。

例子:设 \(G\) 是正方体的对称群,其表现为 \[G = \langle\rho_0,\rho_1,\rho_2\ |\ \rho_0^2=\rho_1^2=\rho_2^2=(\rho_0\rho_1)^4=(\rho_1\rho_2)^3=(\rho_0\rho_2)^2=1\rangle.\] 子群 \(H=\langle \rho_1, \rho_2\rangle\) 是初始顶点的稳定化子群,求 \(G/H\) 的一组右陪集代表元。

我们先罗列一下这个数独游戏已知的关系:

我们的已知关系

  1. \(H\) 的任何生成字 \(w\)\(H\cdot w=H\),即 \(H\rho_1=H\rho_2=H\)。注意此关系仅要求对 \(H\) 成立。
  2. 对任何陪集 \(K\)\(G\) 的任何生成关系 \(r\)\(K\cdot r=K\),即 \(K\rho_i^2=K, i=0,1,2\) 以及 \(K(\rho_0\rho_1)^4=K(\rho_1\rho_2)^3=K(\rho_0\rho_2)^2=K\)。注意此关系要求对所有陪集成立。

这些关系可以存储在两个列表里面,每个关系用一个数组表示。

第一个列表存储的是 \(H\) 的生成字,即

\(H\) 的生成字列表

  1. (1,)
  2. (2,)

第二个列表存储的是 \(G\) 的生成关系,即

\(G\) 的生成关系列表

  1. (0, 0)
  2. (1, 1)
  3. (2, 2)
  4. (0, 1, 0, 1, 0, 1, 0, 1)
  5. (1, 2, 1, 2, 1, 2)
  6. (0, 2, 0, 2)

其中每条关系前面的数字是我们加上的编号以便于称呼。

注:在非 Coxeter 群的情形 (比如后面 snub cube 的情形),还要把每个生成元的逆也作为生成元加入,其在 \(T\) 中也占据一列,所以实际上 \(T\) 的列的个数要 \(\times2\)。但是在 Coxeter 群的情形每个生成元是 2 阶的,其逆元素等于自身,所以不需要额外考虑逆元素。

初始时刻表格 \(T\) 是这样的:

\(\rho_0\) \(\rho_1\) \(\rho_2\)
\(H_0\)

其中 \(H_0\) 代表 \(H\) 对应的陪集。程序首先验证 \(H_0\) 所在的行满足第一个关系列表 (\(H\) 的生成字列表,随后此列表可以被丢弃),然后依次从上到下扫描 \(T\) 的每一行,假设当前扫描的是第 \(i\) 行,对应的陪集为 \(H_i\),程序验证确保对第二个列表 (\(G\) 的生成关系列表) 中的每条关系 \(w\)\(H_i\) 满足 \(H_iw=H_i\),这个过程中可能发现新的陪集,也可能发现已有的某些陪集是重复的,也有可能需要强行定义新的陪集来使得这个验证能够完成。


我们来开始扫描 \(H_0\) 所在的行:首先检查第一个列表中的关系,这个列表仅在扫描 \(H_0\) 时使用一次,扫描完就可以丢弃

(1). 对第 0 条关系 \(H_0\rho_1=H_0\),即 \(T[0][1]=0\)。对第 1 条关系 \(H_0\rho_2=H_0\),即 \(T[0][2]=0\),这时 \(T\) 变成了

\(\rho_0\) \(\rho_1\) \(\rho_2\)
\(H_0\) 0 0

第一个列表扫描完毕,接下来扫描第二个列表。

(2). 对第 2 条关系 \(H_0\rho_0^2=H_0\),由于 \(H_0\rho_0\) 还不知道,我们将其定义为新陪集 \(H_1\) 并将 1 填入 \(T[0][0]\) 位置,此外还要为 \(H_1\) 开辟新的一行:

\(\rho_0\) \(\rho_1\) \(\rho_2\)
\(H_0\) 1 0 0
\(H_1\) 0

注意每次定义新陪集时,比如将 \(H_i\rho_j\) 定义为 \(H_k\),我们同时自动得到了 \(H_k\rho_j=H_i\),因此每次填表时我们都填写一对数字而不是一个

(3). 第 3 条和第 4 条关系已经满足,继续。

(4). 第 5 条关系,\(H_0\rho_0\rho_1\rho_0\rho_1\rho_0\rho_1\rho_0\rho_1=H_0\),我们已经知道 \(H_0\rho_0=H_1\) 但是 \(H_1\rho_1\) 还不知道,将其定义为 \(H_2\),于是 \(T\) 中又添两项,并开辟新的一行给 \(H_2\)

\(\rho_0\) \(\rho_1\) \(\rho_2\)
\(H_0\) 1 0 0
\(H_1\) 0 2
\(H_2\) 1

但是 \(H_2\rho_0\) 还是不知道,所以继续定义 \(H_2\rho_0=H_3\),于是 \(T\) 变成

\(\rho_0\) \(\rho_1\) \(\rho_2\)
\(H_0\) 1 0 0
\(H_1\) 0 2
\(H_2\) 3 1
\(H_3\) 2

但是 \(H_3\rho_1\) 还是不知道,你可能会想把它继续定义为新的陪集 \(H_4\),然后继续扫描。这样做不是不可以,但是这样每次都定义新陪集会生成大量重复的陪集,导致 \(T\) 增长的非常快,对更复杂的群的情形非常耗费计算资源。我们采用更聪明的办法:倒着扫描整个关系,即从右到左扫描 \(H_0\rho_0\rho_1\rho_0\rho_1\rho_0\rho_1\rho_0\rho_1=H_0\) 这条关系。记住我们现在已经正向 (从左到右) 扫描到了下面的位置: \[H_0\rho_0\rho_1\rho_0(=H_3)\rho_1\rho_0\rho_1\rho_0\rho_1=H_0.\] 反向扫描的话很容易看出 \(H_0\rho_1\rho_0\rho_1\rho_0=H_3\),即 \[H_0\rho_0\rho_1\rho_0(=H_3)\rho_1 =H_0\rho_1\rho_0\rho_1\rho_0=H_3.\] 从而 \(H_3\rho_1=H_3\)。这样通过反向扫描我们就推断出了 \(H_3\rho_1\) 的值,避免了定义冗余的陪集。按照 Holt 书中的说法这叫做一个 deduction。这时 \(T\)

\(\rho_0\) \(\rho_1\) \(\rho_2\)
\(H_0\) 1 0 0
\(H_1\) 0 2
\(H_2\) 3 1
\(H_3\) 2 3

所以在实际的程序中,我们总是从关系的两头同时开始扫描,直到它们相遇为止。

(5). 关系 6 已经满足,继续。

(6). 对关系 7 \(H_0\rho_0\rho_2\rho_0\rho_2=H_0\),从两头扫描我们得到 \[H_0\rho_0(=H_1)\rho_2=H_0\rho_2\rho_0=H_1,\]\(H_1\rho_2=H_1\),我们又得到了一个 deduction,从而 \(T\) 变成

\(\rho_0\) \(\rho_1\) \(\rho_2\)
\(H_0\) 1 0 0
\(H_1\) 0 2 1
\(H_2\) 3 1
\(H_3\) 2 3

至此对 \(H_0\) 的扫描全部完成,我们转入扫描 \(H_1\) 所在的行。


注意:从现在起至程序结束,我们不再使用第一个列表

下面开始扫描 \(H_1\) 所在的行。

(1). 经检查关系 2, 3, 4, 5 已经满足,继续。

(2). 对关系 6 有 \(H_1\rho_1\rho_2\rho_1\rho_2\rho_1\rho_2=H_1\),其中 \(H_1\rho_1=H_2\) 已知但 \(H_2\rho_2\) 未知。反向的扫描也会卡在这里: \[H_1\rho_1(=H_2)\rho_2\rho_1=H_1\rho_2\rho_1\rho_2=H_2\rho_2.\] 所以我们定义新陪集 \(H_2\rho_2=H_4\),于是 \(H_4\rho_1=H_4\),从而此时 \(T\)

\(\rho_0\) \(\rho_1\) \(\rho_2\)
\(H_0\) 1 0 0
\(H_1\) 0 2 1
\(H_2\) 3 1 4
\(H_3\) 2 3
\(H_4\) 4 2

(3). 关系 7 已经满足,从而 \(H_1\) 检查完毕,接下来开始扫描 \(H_2\) 的行。


下面开始扫描 \(H_2\) 的行。

(1). 经检查关系 2, 3, 4, 5, 6 都已经满足,继续。

(2). 对关系 7 有 \(H_2\rho_0\rho_2\rho_0\rho_2=H_2\),两边同时扫描的结果为: \[H_2\rho_0(=H_3)\rho_2\rho_0=H_2\rho_2=H_4.\]\(H_3\rho_2\rho_0=H_4\),但是继续正向扫描 \(H_3\rho_2\) 不知道,继续反向扫描 \(H_4\rho_0\) 不知道。定义新陪集 \(H_3\rho_2=H_5\),于是 \(H_5\rho_0=H_4\),我们又可以填入两对 4 个数字,此时 \(T\) 为:

\(\rho_0\) \(\rho_1\) \(\rho_2\)
\(H_0\) 1 0 0
\(H_1\) 0 2 1
\(H_2\) 3 1 4
\(H_3\) 2 3 5
\(H_4\) 5 4 2
\(H_5\) 4 3

\(H_2\) 扫描完毕,下面扫描 \(H_3\) 的行。


我把 \(H_3, H_4, H_5\) 的扫描过程留给你作为练习,\(H_3\) 扫描结束后你得到的 \(T\) 应该如下图:

\(\rho_0\) \(\rho_1\) \(\rho_2\)
\(H_0\) 1 0 0
\(H_1\) 0 2 1
\(H_2\) 3 1 4
\(H_3\) 2 3 5
\(H_4\) 5 4 2
\(H_5\) 4 6 3
\(H_6\) 5 6

\(H_4\) 扫描结束后你得到的 \(T\) 应该如下图:

\(\rho_0\) \(\rho_1\) \(\rho_2\)
\(H_0\) 1 0 0
\(H_1\) 0 2 1
\(H_2\) 3 1 4
\(H_3\) 2 3 5
\(H_4\) 5 4 2
\(H_5\) 4 6 3
\(H_6\) 7 5 6
\(H_7\) 6 7

\(H_5\) 的扫描给不出新的信息。

扫描 \(H_6\) 时,关系 2, 3, 4, 5, 6 都已经满足,由关系 7 \(H_6\rho_0\rho_2\rho_0\rho_2=H_6\) 可得 deduction \(H_7\rho_2=H_7\),于是 \(T\) 可以补全为

\(\rho_0\) \(\rho_1\) \(\rho_2\)
\(H_0\) 1 0 0
\(H_1\) 0 2 1
\(H_2\) 3 1 4
\(H_3\) 2 3 5
\(H_4\) 5 4 2
\(H_5\) 4 6 3
\(H_6\) 7 5 6
\(H_7\) 6 7 7

扫描 \(H_7\) 发现所有关系都已经满足。

至此 \(T\) 的空白位置都已经填满,没有新的陪集可以发现,数独游戏结束,这时得到的 \(T\) 就是 \(G/H\) 的最终乘法表。

由此利用宽度优先搜索不难得出陪集间的关系为: \[\begin{array}{l}H_0 = H_0\cdot e,\\ H_1=H_0\cdot\rho_0,\\H_2=H_1\rho_1=H_0\cdot\rho_0\rho_1,\\H_3=H_2\cdot\rho_0=H_0\cdot\rho_0\rho_1\rho_0,\\ H_4=H_2\cdot\rho_2=H_0\cdot\rho_0\rho_1\rho_2,\\ H_5=H_3\cdot\rho_2=H_0\cdot \rho_0\rho_1\rho_0\rho_2,\\ H_6=H_5\cdot\rho_1=H_0\cdot \rho_0\rho_1\rho_0\rho_2\rho_1,\\ H_7=H_6\cdot\rho_0=H_0\cdot \rho_0\rho_1\rho_0\rho_2\rho_1\rho_0.\end{array}\]

从而其一组陪集代表元可以选为 \[\begin{array}{llll}e&\rho_{0}&\rho_{0}\rho_{1}&\rho_{0}\rho_{1}\rho_{0}\\\rho_{0}\rho_{1}\rho_{2}&\rho_{0}\rho_{1}\rho_{0}\rho_{2}&\rho_{0}\rho_{1}\rho_{0}\rho_{2}\rho_{1}&\rho_{0}\rho_{1}\rho_{0}\rho_{2}\rho_{1}\rho_{0}\end{array} \] 这正是我们前面看到过的。

注意:这个例子看似有点长,但还是一个比较简单的例子,其中并没有出现已知陪集重复的情形 (Holt 的书中称之为 coincidence)。这种情形的处理麻烦一些,因为一旦出现重复的陪集,就有可能 "顺藤摸瓜" 找到更多重复的陪集,再进一步发现更多的重复陪集 ...,这时就要立刻暂停扫描转到处理这个 coincidence:开辟一个栈来存放已知的 coincidence,每次弹出一对,将它们合并,然后把新发现的 coincidence 压入栈中。


总结

最后总结一下整个流程:

计算一个多面体/多胞体 \(P\) 的步骤

  1. [初始化] 找到 \(P\) 对应的 Coxeter-Dynkin 图,并由此确定:
    • \(P\) 的对称群 \(G\) 的一个表现;
    • \(P\) 的初始顶点 \(v_0\) 的坐标;
    • 在空间中摆好反射的镜子,求出每个镜面反射代表的线性变换的矩阵。
  2. [计算顶点] 定出 \(v_0\) 的稳定化子群 \(H\),利用 Todd-Coxeter 算法求出 \(G/H\) 的一组陪集代表元,将它们作用在 \(v_0\) 上得出 \(P\) 的全部顶点。
  3. [计算边] 对每个 \(i\),定出 \(P\) 的一条由反射 \(\rho_i\) 生成的边 \(e\) (如果有的话),求出 \(G/\langle\rho_i\rangle\) 的一组陪集代表元,将它们作用在 \(e\) 上得出 \(P\) 的全部类型为 \(i\) 的边。
  4. [计算面] 对每个 \(i\ne j\),定出 \(P\) 的一个由旋转 \(\rho_i\rho_j\) 生成的面 \(f\) (如果有的话) 及其稳定化子群 \(H\) (不同情况下 \(H\) 可能不同),求出 \(G/H\) 的一组陪集代表元,将它们作用在 \(f\) 上得出 \(P\) 的全部类型为 \((i,j)\) 的面。

这个程序最难的部分在于用 Todd-Coxeter 算法计算 \(P\) 的对称群及其商群的乘法结构,本质上就是求解一个有限 Coxeter 群的所有 word 表示。也有别的方法,比如用 Knuth-Bendix 约化算法 (估计效率会很低),或者用 Brink 和 Howlett 的有限生成 Coxeter 群的极小根理论,但是这算法太太太复杂了。


星状多面体的计算

我还不太清楚所有星状均匀多面体的计算方法,但是至少对三维 Kepler-Poinsot 多面体及其截断,四维的 10 种星状正多胞体及其截断,以及其它部分非凸的多面体,可以采用和上面类似的计算方法,唯一的区别在于需要在生成元 \(\{\rho_0,\rho_1,\rho_2\}\) 之间再加上某些生成关系。

这里以上面的 small-stellated-dodecahedron 为例来说明:其 Coxeter-Dynkin 图为

于是三面镜子的法向量夹角分别为 \(\pi-2\pi/5, \pi/2, \pi-\pi/5\)。如果我们仍然沿用以前的分析,会得出其对称群的表现为

\[G = \langle\rho_0,\rho_1,\rho_2\ |\ \rho_0^2=\rho_1^2=\rho_2^2=(\rho_0\rho_1)^5=(\rho_1\rho_2)^5=(\rho_0\rho_2)^2=1\rangle.\]

这是一个无限群,而且顶点的稳定化子的商群也是无限的,所以还想按以前的方法计算就行不通了。

实际上我们只要在生成关系中再加上一条 \((\rho_0\rho_1\rho_2\rho_1)^3=1\) 即可,即对称群的表现为

\[G = \langle\rho_0,\rho_1,\rho_2\ |\ \rho_0^2=\rho_1^2=\rho_2^2=(\rho_0\rho_1)^5=(\rho_1\rho_2)^5=(\rho_0\rho_2)^2=(\rho_0\rho_1\rho_2\rho_1)^3=1\rangle.\]

这个群用记号 \(\{5, 5\, |\, 3\}\) 表示,这是因为多面体的每个顶点处有 5 个面在此相遇,每个面是一个正五角星,每个面的前三条边全出一个三角形。这个规律最早是 Coxeter 在

Regular skew polyhedra in three and four dimensions, and their topological analogues

中注意到的。

剩下的计算与上面凸的情形完全一样。


Snub cube 的计算

我是不是漏掉了 Snub 多面体的情形没有讲?其实如果你理解了上面的内容,Snub 的情形也是不难理解的。我最好还是用 Snub cube 的例子来说明:

snub cube 和 cube 的区别在于它的对称群只包含旋转,我们已经看到 cube 的对称群 \(G\) 的表现为 \[G = \langle\rho_0,\rho_1,\rho_2\ |\ \rho_0^2=\rho_1^2=\rho_2^2=(\rho_0\rho_1)^4=(\rho_1\rho_2)^3=(\rho_0\rho_2)^2=1\rangle.\] 它有 48 个元素,其中的一半,即 24 个是旋转。这些旋转可以由三个基本旋转 \(r_0=\rho_0\rho_1, r_1=\rho_1\rho_2,r_2=\rho_0\rho_2\) 生成 (由于 \(r_0r_1=r_2\) 因此实际上可以由 \(r_0\)\(r_1\) 生成)。这 24 个旋转就构成了 Snub cube 的对称群 \(\widetilde{G}\)

不难写出 \(\widetilde{G}\) 的表现为 \[\widetilde{G}=\langle r_0,r_1\ |\ r_0^4=r_1^3=(r_0r_1)^2=1\rangle.\]

利用 Todd-Coxeter 算法不难求出这个群的所有 24 个元素:

\[\begin{array}{lll}e&r_{0}&r_{0}r_{0}\\r_{0}r_{0}r_{0}&r_{1}&r_{1}r_{1}\\r_{0}r_{1}&r_{0}r_{1}r_{1}&r_{0}r_{0}r_{1}\\r_{0}r_{0}r_{1}r_{1}&r_{0}r_{0}r_{0}r_{1}&r_{1}r_{0}\\r_{1}r_{0}r_{0}&r_{1}r_{0}r_{0}r_{0}&r_{1}r_{1}r_{0}\\r_{1}r_{1}r_{0}r_{0}&r_{0}r_{1}r_{1}r_{0}&r_{0}r_{1}r_{1}r_{0}r_{0}\\r_{0}r_{0}r_{1}r_{1}r_{0}&r_{1}r_{0}r_{0}r_{1}&r_{1}r_{0}r_{0}r_{1}r_{1}\\r_{1}r_{0}r_{0}r_{0}r_{1}&r_{1}r_{1}r_{0}r_{0}r_{1}&r_{0}r_{1}r_{1}r_{0}r_{0}r_{1}\end{array} \]

将它们作用在任一不属于三个镜面 \(m_0,m_1,m_2\) 的初始点 \(v_0\) 上即得 snub cube 的所有顶点。

snub cube 的边可以这样求:\(r_i(i=0,1,2)\) 作用在 \(v_0\) 上可得一个类型为 \(i\) 的边 \(e_i=(v_0, v_0\cdot r_i)\),将整个 \(G\) 作用在 \(e_i\) 上可得所有类型为 \(i\) 的边。

snub cube 的面可以这样求:由于 \(r_0^4=1\) 所以 \(r_0\) 可以生成一个正四边形的面。类似地由于 \(r_1^3=1\) 所以 \(r_1\) 可以生成一个正三角形的面。还有一种类型的三角形的面可以由 \(r_1\)\(r_2\) 得到。


附录

在项目代码中有一个脚本 run_coset_enumeration.py,可以用来计算一个有限表现群的陪集表 (假定 \(|G/H|<\infty\))。它被设定为读取一个 yaml 文件,该 yaml 文件描述了群 \(G\) 的表现和子群 \(H\) 的生成字,一份样例格式为

name: G8723
relators:
  - aaaaaaaa
  - bbbbbbb
  - abab
  - AbAbAb
subgroup-generators:
  - aa
  - Ab

这里约定大写字母表示对应小写字母的逆元素,如 \(A=a^{-1},B=b^{-1}\) 等等。

于是这个群的表现为 \[G = \langle a, b\ |\ a^8=b^7=(ab)^2=(a^{-1}b)^3=1\rangle.\] 子群 \(H=\langle a^2, a^{-1}b\rangle\)

将此文件存为 filename.yaml 后运行

python run_coset_enumeration.py filename.yaml

输出为

           a    A    b    B
------------------------------
    1:     2    2    3    2
    2:     1    1    1    4
    3:     4    5    6    1
    4:     7    3    2    8
  ...    ...  ...  ...  ...
  ...    ...  ...  ...  ...
  ...    ...  ...  ...  ...
  446:   444  444  441  430
  447:   438  433  432  443
  448:   445  445  440  445

中间我省略了一些输出。于是 \(G/H\) 共有 448 个陪集。

最后附上一张我用此程序制作的 github 头像:(a runcinated 120-cell)


补充

Python 的符号计算包 sympy 里面有一个和这里原理类似的 Todd-Coxeter 的实现,但是那个实现是错误的,对 Mathieu 群 \(M_{12}\) 计算失效,且截至目前 (2018-09-06) 仍未改正,见我提的 issue。我本打算提个 pull request 修正这个错误的,无奈原代码太过冗长啰嗦,多了很多不必要的数据结构和流程,看得人心烦 (实在不懂那个作者的脑回路)。如果做 coset enumeration 计算的话,我个人推荐使用 GAP 及其 ITC 包,那里的实现是这个领域顶尖专家写的,效率和可靠性上都有保障。


2018/11/08 更新

补充一点,可以用一种简洁的途径将一个 3d 多面体渲染在球面上而不必另写专门的代码:直接在 4d 空间中作计算然后投影到 3d 并复用多胞体的渲染代码即可,如下图:


  1. 以本文介绍的知识,这里似乎应该说 \(H\) 保持 \(v_0\) 不变,从而 \(v_0\) 的稳定化子群包含 \(H\),你怎么能断言 \(v_0\) 的稳定化子群就等于 \(H\) 呢?这实际上是 Coxeter 群的一个性质,即在 Coxeter 群 \(W\) 的一个几何实现中 (即把 \(W\) 表示为 \(\mathbb{R}^n\) 中若干关于超平面的反射),对其基本区域闭包中的任何一点 \(v\)\(v\) 的稳定化子群是一个标准椭圆子群,其生成元恰好由超平面包含 \(v\) 的那些单反射组成。这个看似显然的结论证明起来并不平凡,见 Humphreys 的书第一章。

  2. 解释与注解 1 类似。边 \(e\) 的稳定化子群 \(H\) 必然保持 \(e\) 的中点不变,而 \(e\) 的中点必然落在镜面 \(i\) 上,但不会落在其它镜面上,故其稳定化子群就是由单反射 \(\rho_i\) 生成的标准椭圆子群。注意 \(\rho_i\) 虽然保持 \(e\) 不变,但是把有向数组 \((v,u)\) 变为 \((u,v)\),所以要剔除这种重复的边。

  3. 解释与注解 1,2 类似。