Coxeter 群,有限状态机与均匀铺砌

终于完成了一个新的 Python 小程序,虽然还有一些功能有待实现,但是大致效果已经出来了。这个程序是我最近半年多来利用几乎所有业余时间呕心沥血、披星戴月、艰苦攻关完成的新作品,是一个纯粹的、优雅的、有益于广大人民欣赏数学之美的程序,代码在 github 上。

这个程序可以用来绘制二维和三维欧式空间、双曲空间和球面上的各种均匀铺砌 (均匀的含义是对称群在顶点集合上的作用是传递的)。部分例子如下:

例子

  • 下图是二维的欧式铺砌 bitruncated (6, 2, 3),你可以看到每个顶点都有一个编号,编号的规则我在后面会解释。


  • 下图是二维的欧式铺砌 omnitruncated (4, 2, 4):


  • 下图是二维的 Poincaré 圆盘上的正双曲铺砌 (2, 3, 13),绘制了 3447 个顶点,6971 条边,3549 个多边形:


  • 下图是二维的 Poincaré 圆盘上的正双曲铺砌 (3, 8, 2):


  • 下图是二维的 Poincaré 圆盘上的双曲铺砌 omnitruncated (4, 2, 5):


  • 下图是二维的 Poincaré 圆盘上的双曲铺砌 truncated (2, 3, 7):


  • 下图是三维的 Poincaré 单位球中的正双曲铺砌 (3, 5, 3):


  • 下图是三维的 Poincaré 单位球中的正双曲铺砌 (5, 3, 5):


  • 下图是三维的 Poincaré 单位球中的正双曲铺砌 (5, 3, 4):


  • 下图是三维的 Poincaré 单位球中的正双曲铺砌 (4, 3, 5):


  • 以上四个铺砌是三维双曲空间中仅有的正且每个胞腔为紧集的双曲铺砌。如果去掉每个胞腔为紧集的限制的话,但是要求每个胞腔的体积有限的话,则还有 10 个正双曲铺砌,如下图的 (6, 3, 3):

    你可以看到每个胞腔都有一个理想顶点,它位于空间的无穷远处 (即单位球面上)。这 10 个铺砌也叫做仿紧的 (paracompact),因为每个胞腔的体积仍然是有限的。如果再去掉每个胞腔的体积有限的条件的话,得到的铺砌叫做非紧的,这样的正铺砌有无限多个。

  • 如果去掉正铺砌这个限制,只考虑均匀铺砌的话,那例子就非常多了。比如我们可以从正铺砌出发通过截断操作得到许多均匀铺砌,如 rectified (3, 5, 3) 和 rectified (5, 3, 4):

    还有来自 [5, 31,1] 的铺砌 cantic order-5 cubic


  • 下图是一个二维的球面铺砌 omnitruncated (5, 2, 3),绘制在了三维的球面上:


  • 最后是一个 shader 程序,仍然来自 Matt Zucker 的 shadertoy 项目


对二维的铺砌,程序会计算图中每个顶点的坐标、边的下标 (边的下标形如 \((i,j)\)\(i,j\) 是整数,即第 \(i\) 个顶点和第 \(j\) 个顶点构成一条边)、面的下标 (含义类似),然后输出 svg 格式的图形 (这里转成了 png)。对三维的铺砌,程序会计算每个顶点、边、面、胞腔的坐标,然后导出到 POV-Ray 中进行渲染。生成一张二维的铺砌图一般只需几秒或者十几秒,生成一张三维的双曲铺砌则一般需要先花 10 分钟左右计算场景中的 150 万条边 (150 万是我试验发现保证场景中出现足够多的边的一个合适数字),然后花两个小时左右渲染。我的公司办公室有一台 16G 内存的台式机,我一般是下班前运行 Python 的部分来生成场景数据,开启渲染,然后第二天早上上班后查看效果。

图像的绘制原理

网上所有绘制均匀铺砌的程序其基本原理都可以归结为 Wythoff 构造法,这个程序也不例外。即在空间中摆放若干反射平面 (镜子),然后从一个初始点出发,依次关于这些平面作反射变换,得到一些镜像点,再将每个镜像点依次关于这些平面作反射得到新的镜像点,再对新的镜像点重复此步骤,一直下去,就得到了铺砌的全部顶点 (关于 Wythoff 构造法的更多内容可以参考我的另一篇文章)。这个方法在实现时需要维护一个已发现的顶点集合,每次反射后将得到的顶点与集合中元素比较去重 (浮点数向量比较相等一般是取小数点后前若干位进行比较)。这个方法的优点是非常容易实现,可以快速获得顶点和边的坐标 (特别是使用 C++/Java/C# 等编译语言),但是缺点也很明显:

  1. 不容易计算哪些顶点构成一个面?哪些顶点构成一个胞腔?一种解决方案是开始时从一个面/胞腔出发,每次对整个面/胞腔的所有点进行反射,比较去重 (维护一个由面/胞腔的中心组成的集合,中心重合即判定重复),但这是一个很笨的方法。

  2. 这个方法每次遍历每个新的镜像点,将其关于所有平面作反射,实际上得到的结果里面包含了大量的重复点,效率很低。

  3. 浮点数比较去重的过程依赖于机器精度。虽然现在的个人计算机精度已经很高,这样做基本没有什么问题,但从逻辑上讲是一个隐患。想象一下经过几千次反射计算以后由于浮点数误差扩散,导致两个不同的顶点在截取精度以后结果相同呢?

  4. 暴力搜索缺乏数学上的美感。

采用这种方法的一个典型程序见这里

另一种普遍采用的方法叫做 "逆像素计算法",这个方法其实根本不关心顶点坐标和边、面关系之类的任何信息,它只关心最终图像的每个像素的颜色是什么。首先对图像的每个像素点,将其对应到空间中的一点,然后反复进行反射计算,直到反射的结果落入一个事先划定的 "基本区域" (即反射镜面所夹的凸锥) 为止,然后根据最终得到的点在基本区域内的位置对原像素染色。这个方法的优点也是容易实现,可以做出非常精美的图像来。由于它是逐像素染色的,所以适合并行计算,特别是在 shader 里面实现。

这个方法虽然不需要比较去重,但是同样需要作大量的反射计算,而且为了获得高质量的结果需要对图像进行较高的超采样,这更增加了计算量,对 Python 来说是不可承受之重。

采用这种方法的一个典型作品见这里

下面的视频演示的是 Wythoff 构造法的效果,它和万花筒的原理是一样的。我们在二维的 Poincaré 双曲圆盘中的一个房间四周放置若干镜子,则房间中的场景就会在镜子中反复成像,仿佛填满了整个空间。但实际上这个场景里面只有一个初始房间是真实的 (即相机所在的这个),其它的房间都是它在镜子中的虚像。


总之在 Wythoff 构造法中,最关键的部分在于如何快速地、不重复地计算所有顶点坐标,以及它们之间的边、面、胞腔的关系。这正是这个程序的优雅之处,它采用的是群论的方法,首先在铺砌的对称群中进行符号计算,得出每个顶点在最短字典序下对应的单词表示,边的下标,面的下标等信息,然后将每个顶点对应的单词作用在初始顶点上得到该顶点的坐标。换句话说,在计算每个顶点的最终坐标之前,它已经提前计算好了有多少个顶点,每个顶点是由初始顶点经过哪些反射得到的,哪些顶点构成边,哪些顶点构成面,哪些顶点构成胞腔,等等。这些计算仅涉及整数运算,完全避免了浮点数精度损失的问题。是不是很神奇呢?

下面我用一个具体的例子来演示这个算法的步骤。

例子:omnitruncated (7, 2, 3) 铺砌

Omnitruncated (7, 2, 3) 铺砌对应的 Coxeter-Dynkin 图如下:

这是一个二维的双曲铺砌 1,其对称群是由 Coxeter 矩阵

\[M=\begin{pmatrix} 1 & 7 & 2 \\ 7 &1 &3\\ 2 & 3 &1\end{pmatrix}\]

确定的 Coxeter 群

\[G = \langle s_0,s_1, s_2\ |\ s_0^2=s_1^2=s_2^2=(s_0s_1)^7=(s_1s_2)^3=(s_0s_2)^2=1\rangle.\]

初始顶点不包含于任何镜面上,所以其稳定化子群只包含单位元 1,从而根据轨道—稳定化子定理 \(G\) 的每一个元素都对应铺砌中的一个顶点。

\(G\) 中每个元素 \(g\) 都可以表示为生成元 \(s_0,s_1,s_2\) 的乘积,我们称任何这样的表示方式为 \(g\) 的一个单词表示。\(g\) 的单词表示不是唯一的,比如从定义 \(G\) 的生成关系中可以看出 \(s_0s_2=s_2s_0\), \(s_1s_2s_1=s_2s_1s_2\) 等等。但是在 \(g\) 的所有单词表示中,我们总是可以选择 "最小的" 那个来作为 \(g\)规范表示,这个排序的依据就是最短字典序 (shortlex order):首先规定生成元 \(s_0,s_1,s_2\) 之间的字母顺序关系为 \(s_0<s_1<s_2\),然后将这个顺序扩展到任何两个单词 \(w_1\)\(w_2\) 上去:

定义:设 \(w_1=s_{i_1}s_{i_2}\cdots s_{i_n}\)\(w_2=s_{j_1}s_{j_2}\cdots s_{j_m}\) 是两个不同的单词,它们二者之间的顺序关系由如下方式确定:

  1. 先比较长度,长度长的为较大者,即若 \(n<m\)\(w_1<w_2\)\(n>m\)\(w_1>w_2\)
  2. 若长度相等则按字母顺序从左到右逐项比较,设 \(k\) 使得对任何 \(l<k\)\(s_{i_l}=s_{j_l}\)\(s_{i_k}\ne s_{j_k}\),则规定 \(w_1,w_2\) 之间的大小关系与 \(s_{i_k},s_{j_k}\) 之间的大小关系相同。

由于最短字典序是所有由 \(s_0,s_1,s_2\) 组成的单词集合上的全序,在这个全序下任何 \(g\in G\) 都有唯一的规范表示,任何两个 \(g,g'\in G\) 之间的大小顺序就规定为它们的规范表示之间的大小关系。一个元素 \(g\in G\) 的长度定义为其规范表示的长度。

定义 \(\mathcal{SL}(G)\)\(G\) 中所有元素的规范表示组成的集合,下面列出了 \(\mathcal{SL}(G)\) 中所有长度不超过 5 的元素,一共有 37 个:(从小到大按行排列)

\[ \begin{array}{lllll}e&s_{0}&s_{1}&s_{2}&s_{0}s_{1}\\s_{0}s_{2}&s_{1}s_{0}&s_{1}s_{2}&s_{2}s_{1}&s_{0}s_{1}s_{0}\\s_{0}s_{1}s_{2}&s_{0}s_{2}s_{1}&s_{1}s_{0}s_{1}&s_{1}s_{0}s_{2}&s_{1}s_{2}s_{1}\\s_{2}s_{1}s_{0}&s_{0}s_{1}s_{0}s_{1}&s_{0}s_{1}s_{0}s_{2}&s_{0}s_{1}s_{2}s_{1}&s_{0}s_{2}s_{1}s_{0}\\s_{1}s_{0}s_{1}s_{0}&s_{1}s_{0}s_{1}s_{2}&s_{1}s_{0}s_{2}s_{1}&s_{1}s_{2}s_{1}s_{0}&s_{2}s_{1}s_{0}s_{1}\\s_{0}s_{1}s_{0}s_{1}s_{0}&s_{0}s_{1}s_{0}s_{1}s_{2}&s_{0}s_{1}s_{0}s_{2}s_{1}&s_{0}s_{1}s_{2}s_{1}s_{0}&s_{0}s_{2}s_{1}s_{0}s_{1}\\s_{1}s_{0}s_{1}s_{0}s_{1}&s_{1}s_{0}s_{1}s_{0}s_{2}&s_{1}s_{0}s_{1}s_{2}s_{1}&s_{1}s_{0}s_{2}s_{1}s_{0}&s_{1}s_{2}s_{1}s_{0}s_{1}\\s_{2}s_{1}s_{0}s_{1}s_{0}&s_{2}s_{1}s_{0}s_{1}s_{2}&\end{array} \]

注意由 \(s_0,s_1,s_2\) 组成的所有长度不超过 5 的单词 (含单位元) 有 \(1+3+\cdots+3^5=364\) 个,上表告诉我们实际上它们只包含了 \(G\) 中 37 个不同的元素,即将它们作用在铺砌中的初始顶点上只会得到 37 个顶点,其余 364 - 37 = 327 个都是重复的。进一步计算可得长度不超过 6 的单词有 1093 个而实际上只包含了 53 个不同的元素。所以如果我们能够快速地生成 \(\mathcal{SL}(G)\) 中的元素而不是去遍历所有可能单词的话,就可以大大提高计算效率。

那么怎么才能"快速精准地" 生成 \(\mathcal{SL}(G)\) 中的元素呢?这就引出了一个关于 Coxeter 群的重要结论:

定理 [Brigitte Brink & Robert B. Howlett, 1993]:若 \(G\) 是有限生成 Coxeter 群,则 \(\mathcal{SL}(G)\) 是一个正则语言

正则语言这个术语来自计算机科学,关于正则语言的一个基本事实是,一个有限字符集上的正则语言总是可以被一个确定的有限状态自动机 (definite finite automaton) 识别,这样的有限状态机不是唯一的,但是在等价意义下总是存在一个最小的 (含有的状态数目最少),下图显示的是识别 (7, 2, 3) 这个群的 \(\mathcal{SL}(G)\) 的最小状态机:(下图在 graphviz 输出的基础上进行了二次绘制)

你可以看到在上图中一共有 19 个节点 (即状态),每个状态都有一个编号,这个编号没有任何意义,可以不用理会,实际上给状态重新编号不影响有限状态机识别的语言。圈红的节点 0 是初始状态。

图中的每条有向边规定了状态之间的转移规则,边的标号 \(i\) 代表生成元 \(s_i\)。从初始状态出发,每次沿着一条有向边移动到下一个状态,经过的路径给出了一个由 \(s_0,s_1,s_2\) 组成的单词,所有的路径给出的单词组成的集合就是这个有限状态机识别的语言,即 \(\mathcal{SL}(G)\)

举个例子:

  1. 长度为 0 的路径对应的是 \(G\) 的单位元。
  2. 长度为 1 的三条路径 \(0\xrightarrow{\ s_0\ }1\), \(0\xrightarrow{\ s_1\ }2\), \(0\xrightarrow{\ s_2\ }8\) 对应的是 \(G\) 的三个长度为 1 的元素 \(s_0,s_1,s_2\)
  3. 长度为 2 的 5 条路径 \(0\xrightarrow{\ s_0\ }1\xrightarrow{\ s_1\ }2\), \(0\xrightarrow{\ s_0\ }1\xrightarrow{\ s_2\ }8\), \(0\xrightarrow{\ s_1\ }2\xrightarrow{\ s_0\ }3\), \(0\xrightarrow{\ s_1\ }2\xrightarrow{\ s_2\ }8\), \(0\xrightarrow{\ s_2\ }8\xrightarrow{\ s_1\ }9\) 对应的是 \(G\) 的 5 个长度为 2 的元素 \(s_0s_1,s_0s_2,s_1s_0,s_1s_2,s_2s_1\)

以此类推,我们可以很容易用宽度优先搜索来不重复不遗漏地遍历任意长度范围内的群元素。特别地如果你 "按图索骥" 验证一下的话可以发现所有长度不超过 5 的路径一共有 37 条,它们正对应前面列出的 \(\mathcal{SL}(G)\) 中长度不超过 5 的 37 个单词。

注意无限 Coxeter 群的有限状态机必然含有回路,而有限 Coxeter 群的有限状态机则必然是一个以初始状态为根节点的有向树,例如下图显示的是正四面体的对称群 \(S_4\) 的有限状态机:


不难遍历得出其 24 个不同的路径,它们分别对应 \(S_4\) 的 24 个元素的最短字典序表示:

\[ \begin{array}{llll}e&s_{0}&s_{1}&s_{2}\\s_{0}s_{1}&s_{0}s_{2}&s_{1}s_{0}&s_{1}s_{2}\\s_{2}s_{1}&s_{0}s_{1}s_{0}&s_{0}s_{1}s_{2}&s_{0}s_{2}s_{1}\\s_{1}s_{0}s_{2}&s_{1}s_{2}s_{1}&s_{2}s_{1}s_{0}&s_{0}s_{1}s_{0}s_{2}\\s_{0}s_{1}s_{2}s_{1}&s_{0}s_{2}s_{1}s_{0}&s_{1}s_{0}s_{2}s_{1}&s_{1}s_{2}s_{1}s_{0}\\s_{0}s_{1}s_{0}s_{2}s_{1}&s_{0}s_{1}s_{2}s_{1}s_{0}&s_{1}s_{0}s_{2}s_{1}s_{0}&s_{0}s_{1}s_{0}s_{2}s_{1}s_{0}\end{array} \]

问题 1:怎样计算 \(\mathcal{SL}(G)\) 对应的有限状态机?

此问题技术太过复杂,完整介绍全部内容的话本文难以承受,后面附有一个简单的解释。我当时主要参考了 Casselman 的讲义 2 3 4,以及 Humphreys 的教材 5。我认为这几份文献的内容足够了。

有了每个群元素的唯一的规范表示,我们就可以很容易地计算铺砌中每个顶点的坐标了:

\(w=s_{i_0}s_{i_1}\cdots s_{i_n}\),初始顶点为 \(v_0\),则 \(w\)\(v_0\) 上的作用为 \[w\cdot v_0 = s_{i_0}(s_{i_1}(\cdots s_{i_n}(v_0)).\] 即从右到左依次计算每个生成元作用的结果。当然由于 \(G\) 是无限群,我们只能计算那些长度不超过一定范围的群元素对应的顶点。假设我们已经有了前面这 37 个顶点,把它们按照最短字典序排序存储在一个列表 \(L\) 里。为了绘制铺砌中的边,我们还需要计算 \(L\) 中哪些顶点是相邻的,这个怎么解决呢?

首先我们需要计算一个 \(L\) 中元素之间的乘法表 \(T\)\(T\) 的作用是帮助我们查找任何一个单词,其规范表示在 \(L\) 中的下标 (不在 \(L\) 中的话则返回 None)。\(T\) 是一个大小为 \(|L|\times 3\) 的二维数组,其第 \(i\) 行对应 \(L\) 中的第 \(i\) 个群元素 \(w_i\),第 \(j\) 列对应群的生成元 \(s_j\)\(T[i][j]\) 的值等于 \(s_jw_i\) 这个元素 (注意这个乘积未必是个规范表示) 的规范表示在 \(L\) 中的下标。如果这个元素不在 \(L\) 中则以 None 代替。

在我们的例子中 \(T\) 的值如下,\(L\) 中的元素放在了第二列:

点击展开列表 \(T\)
V word \(s_0\) \(s_1\) \(s_2\)
0 \(e\) 1 2 3
1 \(s_{0}\) 0 6 5
2 \(s_{1}\) 4 0 8
3 \(s_{2}\) 5 7 0
4 \(s_{0}s_{1}\) 2 12 11
5 \(s_{0}s_{2}\) 3 13 1
6 \(s_{1}s_{0}\) 9 1 15
7 \(s_{1}s_{2}\) 10 3 14
8 \(s_{2}s_{1}\) 11 14 2
9 \(s_{0}s_{1}s_{0}\) 6 20 19
10 \(s_{0}s_{1}s_{2}\) 7 21 18
11 \(s_{0}s_{2}s_{1}\) 8 22 4
12 \(s_{1}s_{0}s_{1}\) 16 4 24
13 \(s_{1}s_{0}s_{2}\) 17 5 23
14 \(s_{1}s_{2}s_{1}\) 18 8 7
15 \(s_{2}s_{1}s_{0}\) 19 23 6
16 \(s_{0}s_{1}s_{0}s_{1}\) 12 30 29
17 \(s_{0}s_{1}s_{0}s_{2}\) 13 31 28
18 \(s_{0}s_{1}s_{2}s_{1}\) 14 32 10
19 \(s_{0}s_{2}s_{1}s_{0}\) 15 33 9
20 \(s_{1}s_{0}s_{1}s_{0}\) 25 9 35
21 \(s_{1}s_{0}s_{1}s_{2}\) 26 10 36
22 \(s_{1}s_{0}s_{2}s_{1}\) 27 11 34
23 \(s_{1}s_{2}s_{1}s_{0}\) 28 15 13
24 \(s_{2}s_{1}s_{0}s_{1}\) 29 34 12
25 \(s_{0}s_{1}s_{0}s_{1}s_{0}\) 20 None None
26 \(s_{0}s_{1}s_{0}s_{1}s_{2}\) 21 None None
27 \(s_{0}s_{1}s_{0}s_{2}s_{1}\) 22 None None
28 \(s_{0}s_{1}s_{2}s_{1}s_{0}\) 23 None 17
29 \(s_{0}s_{2}s_{1}s_{0}s_{1}\) 24 None 16
30 \(s_{1}s_{0}s_{1}s_{0}s_{1}\) None 16 None
31 \(s_{1}s_{0}s_{1}s_{0}s_{2}\) None 17 None
32 \(s_{1}s_{0}s_{1}s_{2}s_{1}\) None 18 None
33 \(s_{1}s_{0}s_{2}s_{1}s_{0}\) None 19 None
34 \(s_{1}s_{2}s_{1}s_{0}s_{1}\) None 24 22
35 \(s_{2}s_{1}s_{0}s_{1}s_{0}\) None None 20
36 \(s_{2}s_{1}s_{0}s_{1}s_{2}\) None None 21

于是对任意的单词 \(w=s_{i_0}s_{i_1}\cdots s_{i_n}\),我们可以从 \(T\) 的第 0 行出发,先找到 \(s_{i_n}\)\(L\) 中对应的元素,假设是第 \(k\) 个,那么就跳到第 \(k\) 行,由 \(s_{i_{n-1}}\) 对应的列找到 \(s_{i_{n-1}}s_{i_n}\)\(L\) 中对应的元素,再顺藤摸瓜找到 \(s_{i_{n-2}}s_{i_{n-1}}s_{i_n}\)\(L\) 中对应的元素,...,如此下去即可确定 \(w\)\(L\) 中对应的元素 (或者 None)。

假设初始顶点 \(v_0\) 关于第 \(i\) 面镜子反射后得到的虚像是 \(v_1=s_i(v_0)\),则 \(e=(v_0,v_1)\) 构成一条类型为 \(i\) 的边。边 \(e\) 的稳定化子群就是标准椭圆子群 \(H=\langle s_i\rangle\)。根据轨道稳定化子定理,铺砌中所有类型为 \(i\) 的边都可以通过将 \(G/H\) 中某个代表元作用在 \(e\) 上得到。显然 \(e\) 的两个端点对应的规范表示分别是单位元 1 和 \(s_i\),则对任一 \(w\in G\),我们首先计算 \(w\) 关于 \(H\) 的陪集代表元 \(w_H\) (这一步需要配合去重,去掉重复的陪集代表元),\(w_H\cdot e\) 的两个端点对应的单词分别是 \(w_H\)\(w_Hs_i\),然后按照上面的步骤找到它俩在 \(L\) 中对应的元素下标 (如果有一个是 None 说明这条边连接了不在 \(L\) 中的顶点),这就得到了边的下标。

\(L\) 中 37 个顶点构成的边如下:

图中标号 0 的顶点是初始顶点,其对应的单词是单位元 1。边中白色的线条个数表示这条边对应的生成元:0 条表示 \(s_0\),1 条表示 \(s_1\),2 条表示 \(s_2\)

从这个图里可以看出很多有用的信息,非常有助于理解 Coxeter 群的最短字典序表示与铺砌顶点的对应关系:

  1. 首先可以看出将顶点 0 关于三个镜面进行一次反射得到的新虚像是 1, 2, 3,它们对应的单词是 \(L\) 中长度为 1 的 \(s_0,s_1,s_2\)
  2. 将顶点 0 关于各个镜面进行两次反射得到的新虚像是 4, 5, 6, 7, 8,它们对应的单词是 \(L\) 中长度为 2 的 \(s_0s_1,s_0s_2,s_1s_0,s_1s_2,s_2s_1\)
  3. 所有顶点与初始顶点之间的最短路径的长度都不超过 5。
  4. 从图中我们可以看出每个顶点对应的单词的规范表示。例如从 0 号顶点到 33 号顶点有两条最短路径:\[ \begin{align*}&0\xrightarrow{\ s_1\ }2\xrightarrow{\ s_0\ }6\xrightarrow{\ s_2\ }13\xrightarrow{\ s_1\ }22\xrightarrow{\ s_0\ }33.\\ &0\xrightarrow{\ s_1\ }2\xrightarrow{\ s_2\ }7\xrightarrow{\ s_0\ }13\xrightarrow{\ s_1\ }22\xrightarrow{\ s_0\ }33. \end{align*}\] 只要按顺序连起来 (从左到右) 就可以得到 33 号顶点对应的两个单词:\(s_1s_0s_2s_1s_0\)\(s_1s_2s_0s_1s_0\),它们都把 0 号顶点变为 33 号顶点,但是前者才是最短字典序表示。这一点也可以从顶点标号中看出来:从二号顶点开始两条路径分别去了 6 号和 7 号顶点,由于 6 号顶点在 \(\mathcal{SL}(G)\) 中更小,因此这条路径必然才是 33 号顶点规范表示对应的路径。
  5. 一个单词 \(w=s_{i_0}s_{i_1}\cdots s_{i_n}\) 作用在 0 号顶点上 (\(s_{i_n}\) 先作用,\(s_{i_0}\) 在最后) 的结果是一条从 0 号顶点出发,标号依次为 \(s_{i_0}, s_{i_1},\ldots,s_{i_n}\) 的路径 (\(s_{i_0}\) 在先,\(s_{i_n}\) 在最后),与作用的顺序相反。

你可能要问,直接计算 \(w_H\)\(w_Hs_i\) 的规范表示,然后去 \(L\) 中查找下标不就可以了吗?为什么计算这个表 \(T\) 呢?这里我的考虑有两个:

  1. \(T\) 的计算是一次性的,而且查表速度很快。不借助 \(T\) 的话就要每条边都计算两次规范表示然后进行两次 \(L\) 中的逐项比较查找。

  2. \(T\) 也同样用于计算面和胞腔。不借助 \(T\) 的话就要对每个胞腔的每个顶点 (一个胞腔是一个多面体,可能有几十个顶点) 计算一次规范表示并逐项比较查找。

问题 2:给定任一单词 \(w\) 和标准椭圆子群 \(H\),怎样计算 \(w\) 的规范表示?怎样计算 \(w\) 关于 \(H\) 的陪集代表元?

简单的介绍见后,详细证明见 Casselman 的文章,本文暂且按下不表。

计算面的步骤也是完全类似的,初始顶点关于 \(i, j\) 两面镜子的反射的复合是一个旋转,这个旋转连续作用 \(m\) 次可以生成一个正多边形的面,其中 \(m\) 是 Coxeter 矩阵中的 \((i,j)\) 分量。这个多边形的稳定化子群是标准椭圆子群 \(H=\langle i, j\rangle\),我们仍然可以得出每个顶点对应的一个单词表示,用 \(G/H\) 里的代表元作用在上面,然后去 \(L\) 里面查找对应的规范表示的下标。

最终得到的图像如下图,计算了 30517 个顶点,42057 条边,11541 个多边形:

关于代码

代码中最核心的部分是一个叫 coxeter 的模块,它主要实现的是 CoxeterGroup 类,这个类由一个 Coxeter 矩阵初始化,负责计算该 Coxeter 群的极小根反射表、实现群元素的乘法、陪集的计算,计算并绘制有限状态机。在另一个 tiling.py 文件中调用此类,计算铺砌的顶点、边、面并绘制。

绘制 \(\mathcal{SL}(G)\) 的有限状态机需要使用 pygraphviz 模块,这个模块依赖于 graphviz 软件和 libgraphviz-dev

我把铺砌的绘图大部分 "外包" 出去了,二维双曲空间的绘图用的是一个叫 hyperbolic 的三方库,这个库的代码写的挺糟糕的,然而我一时也没有精力分出来再写一个,所以暂且先凑合着。二维情形的绘图主要麻烦的地方在为了绘制具有常 "双曲宽度" 的边需要计算对应的 hypercycle。 球面和三维双曲情形的绘图用的是 POV-Ray,主要使用了 sphere_sweep 这个对象。

算法运行效率

这个算法比较耗时的地方有两个:

  1. 计算群 \(G\)极小根部分。计算 Coxeter 群的极小根是计算 \(\mathcal{SL}(G)\) 的有限状态机、计算群的规范表示所必须的重要步骤,而极小根的计算量依赖于 Coxeter 群的 "复杂度":对 \((p, q, r)\) 这样的三角群,当 \(p,q,r\) 都不大于 10 时计算还比较快 (双曲的话稍慢些,要花个几秒钟),然而对 (19, 20, 21) 这样的群计算速度就很慢了。这一点与逆像素反射方法是不同的,逆像素反射法的计算量几乎不随群的变化而变化。

  2. 计算乘法表 \(T\)。对于三维双曲铺砌要生成一个效果较好的场景至少需要几十万个顶点,对这么多顶点计算乘法表非常耗时,内存和查找开销都很大。这种情形我采用的是有限状态机生成群元素 + 传统的浮点数去重,一边计算一边输出以节省空间开销。

总结起来就是对不太复杂的群,或者需要生成的数据量不太大的时候这种方法的速度还不错。

下一步计划

这篇文章遗漏了许多技术性的细节没有解释,比如怎样计算识别一个 Coxeter 群的最短字典序语言的有限状态机,怎样计算一个单词的规范形式,怎样计算一个单词对应元素关于一个标准椭圆子群的的陪集代表元,还有 Coxeter 群的分类以及对应 Tit 锥与铺砌的关系等等。它们对理解这个程序的代码都是至关重要的。下一步我会写一个系列,详细介绍其中的数学。

程序方面已知需要做如下改进:

  1. 目前的代码比较粗糙,需要重构改进结构和可读性。
  2. 替换掉 hyperbolic 这个第三方库。
  3. 增加酷炫的 shader 动画演示。
  4. 增加绘制具有超理想 (hyperideal) 顶点的三维双曲铺砌。现在比较纠结的是 Python 逐个像素计算这种图像非常费劲,在考虑是否用 C++ 或者 C# 实现。

附录:对若干关键点的解释

在这里我简单介绍计算 \(\mathcal{SL}(G)\) 的有限状态机、计算任一单词的规范形式、计算陪集代表元的步骤。需要你了解 Coxeter 群的几何实现、Tit's cone、Coxeter 群的根系等预备知识,这部分内容可以参考 Humphreys 的教材。

以上所有的计算都基于一个叫做极小根反射表的东西。仍然以 (7, 2, 3) 这个群为例子,先看下图:


这个图和上图一样,只不过多了 12 个标记出的镜面,这 12 个镜面有特殊的含义:它们是群 \(G\) 的根系中的 12 个极小根。

你可以把 \(G\) 的根系理解为图中所有白色的镜面,每个镜面是一条测地线,它交单位圆的边界于两点。这些镜面都是三个初始镜面在群 \(G\) 作用下的结果。每个镜面把单位圆 (即整个双曲空间) 分割成两部分,这两部分分别对应根系中的一个正根和一个负根,其中基本区域三角形 \(\Delta ABC\) 所在的半空间是正根对应的部分。

从直观上说,极小根 \(r\) 是那些满足如下条件的镜面:不存在任何镜面 \(r'\ne r\),使得从 \(\Delta ABC\) 往外看,视线会被 \(r'\) 完全挡住而看不到 \(r\) 的任何部分。单根必然都是极小根 (它们是基本区域的边界)。关于极小根的一个重要事实是,任何有限生成 Coxeter 群的极小根的个数都是有限的。这个结论在 Brink 和 Howlett 的证明中起到最关键的作用。

极小根的反射表 reftable 定义如下:这也是一个二维数组,其第 \(i\) 行对应第 \(i\) 个极小根 \(\alpha_i\),其第 \(j\) 列对应第 \(j\) 个生成元 \(s_j\),其 \((i, j)\) 位置填入的是 \(\beta=s_j(\alpha_i)\) 的结果:

  1. 如果 \(\beta=\alpha_k\) 是第 \(k\) 个极小根,则填入 \(k\)
  2. 如果 \(\beta\) 是一个负根 (此情形发生当且仅当 \(\alpha_i\) 是第 \(j\) 个单根),则填入 -1
  3. 如果 \(\beta\) 是一个正根,但不是极小根,则填入 None

(7, 2, 3) 群的极小根反射表如下:

root \(s_0\) \(s_1\) \(s_2\)
0 -1 3 0
1 4 -1 5
2 2 5 -1
3 6 0 7
4 1 8 9
5 9 2 1
6 3 10 11
7 11 7 3
8 10 4 None
9 5 None 4
10 8 6 None
11 7 None 6

\(G\) 的所有极小根之集为 \(\Sigma\)\(\mathcal{G}\) 的有限状态机的状态都是 \(\Sigma\) 的子集,状态之间的转移关系规定如下:

\[S\xrightarrow{\ s_i\ } \{s_i\} \cup (s_i(S)\cup\{ s_i(\alpha_j),j<i\})\cap\Sigma.\]

由此不难利用宽度优先搜索获得所有的状态以及它们之间的转移关系,然后用 Hopcroft 算法将得到的有限状态机最小化。

下图显示的是 (7, 2, 3) 这个群的状态机中,每个状态对应的 \(\Sigma\) 的子集:

计算两个群元素乘法的规范表示的代码如下,其中 s 是一个生成元,word 是一个逆最短字典序下的单词 (逆最短字典序就是把最短字典序中的单词反过来),函数返回的也是对应结果在逆最短字典序下的规范形式:

def left_mul_invshortlex(reftable, s, word):
word = tuple(word)
t = s
k = -1
mu = s
for i, s_i in enumerate(word):
if mu is None:
return word[:k+1] + (t,) + word[k+1:]
elif mu < 0:
return word[:i] + word[i+1:]
elif mu < s_i:
t = mu
k = i
else:
continue
return word[:k+1] + (t,) + word[k+1:]

由此函数不难计算任何两个单词相乘的逆最短字典序下的规范形式,倒过来自然也就解决了最短字典序下的乘法问题。

这里先在逆字典序下计算然后再倒过来获得字典序是为了和 Casselman 的文章中的论述保持一致。

计算一个规范表示 word 关于某个标准椭圆子群的陪集是最简单的:设 \(H\) 是一个标准椭圆子群,其生成元是 \(G\) 的生成元的一个子集 \(T\),则计算 word 关于 \(H\) 的左陪集代表元的伪代码如下:

begin
x := word
u := 1
while l(xt) < l(x) for some t in T
x := xt
u := tu
end

其中 \(l(\cdot)\) 是 Coxeter 群上的长度函数。这个算法会把 word 分解为形如 \(x^T\cdot w_T\) 的形式,其中 \(w_T\in H\),且对任何 \(t\in T\)\(l(x^Tt)>l(x^T)\)。最终得到的陪集代表元 \(x^T\) 是规范表示。

对有限 Coxeter 群,其所有正根都是极小根;对仿射 Coxeter 群,其根系由一些平行的镜面族构成,每族镜面中的反射镜面互相平行。每个族中存在一对极小根,它们把基本区域 \(\Delta ABC\) 夹在中间并完全挡住本族中外面的镜面,所以只有它俩才是极小根。下图显示的是 (6, 2, 3) (仿射 \(\widetilde{G}_2\)) 的 12 个极小根:


  1. 一个铺砌的类型可以由其 Cartan 矩阵 \(C=-\cos(\pi/M)\) 确定,是球面铺砌当且仅当 \(C\) 是正定的,是欧式铺砌当且仅当 \(C\) 是半正定的,是双曲铺砌当且仅当 \(C\) 是不定的,且 Coxeter-Dynkin 图的任何连通子图给出的铺砌都是球面的或者欧式的。

  2. Automata to perform basic calculations in Coxeter groups, by Bill Casselman.

  3. Computation in Coxeter groups I. Multiplication, by Bill Casselman.

  4. Computation in Coxeter groups II. Constructing minimal roots, by Bill Casselman.

  5. Reflection Groups and Coxeter Groups, by James E. Humphreys.

当前网速较慢或者你使用的浏览器不支持博客特定功能,请尝试刷新或换用Chrome、Firefox等现代浏览器