引言
本文主要关注数独 (Sudoku) 的计数问题，也即合法的数独排列的总数，和考虑到对称性之后，真正独立的解的个数。事实上，数独对应的数学结构，解析性并不好，和数独相关的大部分问题也是 NP 困难的，没有好的解析解答和快捷方法，只能暴力求解。但好在数独的问题尺寸不大，多数时候暴力求解也不很复杂。n=3 （也即常见的 9*9 的数独）的解的个数问题，在 2005 年已经被解决了1。其计算方式的主体也是 brute force 列举。但本文更多的只是想展示这类问题，特别是考虑对称性后独立数独排布计数的方案，因此全文会采取 n=2 的数独进行讲解，有效地降低问题的暴力程度的同时，保留了需要解决 n=3 问题的全部要素。
数独个数
这里我们直接给出我们要考虑的 n=2 的数独的定义，我们需要用 1，2，3，4 填满 4*4 的方格，使得每行，每列，和4个 2*2 的子块，都由四个不同的数字组成。我们想解答的问题是，满足以上限制的数字填法究竟有多少种，其中又究竟有多少种是本质上不一样的。我们先看第一个问题，最简单的想法就是暴力穷举。所有在 4*4 格子中填入 4 个 1，2，3，4 的方案一共有 \(C_{16}^4C_{12}^4C_{8}^4=63063000\) 种，这不到一亿种可能性，每种验证下合法性，记个数，也就是一眨眼的事情。但是本文的 n=2 数独的记数是一个 n=3 数独计数的教学，在 n=3 时，总的排布可能性高达 \(10^{70}\) 的数量级，是无论如何也穷举不来的。
因此我们先考虑置换对称性，由于 1234 的角色任意置换，合法数独还是合法的，因此我们可以只考虑第一个子块是按 1，2，3，4 排列的数独的个数。将此数目乘以 \(4!=24\) 就是总的合法数独的个数。在此基础上，我们在考虑前两行剩余部分，也即第二个子块的数字填充。考虑到我们还有第三列和第四列交换，合法数独依旧合法的对称性，我们可以总是让第一行第三列的数小于第一行第四列的数。也即我们可以实现 1234，34xx，xxxx，xxxx 的构型，并且只需考虑 x 处的合法填充即可。所得的数目再乘以 48 就是总的数独数。这样我们有两种独立的前两行的填充可能，也即第二行的两个 x 添入 12 或者是 21。我们只需分别算出两者对应的下两行的填充可能性即可。这可以通过一个简单的深度优先搜索实现穷举后两行的合法构型。
def possible_value(i, j, st):
pos = set([1, 2, 3, 4])
for n in st[i, :]:
if n != 0:
pos.discard(n)
for n in st[:, j]:
if n != 0:
pos.discard(n)
for n in st[
(i // 2) * 2 : (i // 2) * 2 + 2, (j // 2) * 2 : (j // 2) * 2 + 2
].reshape([4]):
if n != 0:
pos.discard(n)
return pos
def next_id(i, j):
if j < 3:
return i, j + 1
elif j == 3 and i < 3:
return i + 1, 0
else:
raise Exception("illegal index")
def find_solution(i, j, st, r):
if st[i, j] != 0:
return find_solution(*next_id(i, j), st, r)
ps = possible_value(i, j, st)
if i == 3 and j == 3:
for p in ps:
st[i, j] = p
r.append(st.copy())
st[i, j] = 0
return r
for p in ps:
st[i, j] = p
find_solution(*next_id(i, j), st, r)
st[i, j] = 0
return r
r = []
st = np.array([[1,2,3,4],[3,4,0,0],[0,0,0,0],[0,0,0,0]])
sols = find_solution(0, 0, st, r)
print(len(sols)) #6
由此可知 n=2 的数独解的总数为 \(6*2*4!=288\)。n = 3 的数独计数精神上和这边是一样的。首先是把左上角的子块通过置换对称性固定成 1-9。然后考虑前三行另两个子块的填充情况，可以分析得到，前三行总共有 \(2*(3!)^6+18*3*(3!)^6 = 2612736\) 种情况。对于每种情况暴力穷举下面六行的可能情况的话，工作量还是过于大。因此依旧要充分考虑利用保持数独合法性的对称性，比如我们在 n=2 时考虑的列交换等等，只不过 n=3 时需要利用的对称性还有很多。通过这些对称性，我们可以把这两百多万种情况合并成较少的等价类。这里等价类的意思是，每类里不同的前三行填充，对应的后六行合法填充的数目一定是相同的（对称性保证）。通过大量的枚举和对称性分析，可以最后约化到 44 个等价类。其中有些对称性约化比较复杂，但实际上我们不做也行，只不过需要穷举的等价类略多而已，不影响结果的正确性。之后我们只需对于每一个前三行排布的等价类，对后六行做回溯来填充，并计算可能的数目最后求和即可。每种等价类构型对应的前三行排布的数目，和后六行可能的填充情况数可以参考表格2，详细的前三行的对称性规约可以参考3。n = 3 数独总数的结果为 6670903752021072936960，大小为 n 的数独的总数数列为 A107739。
独立数独个数
本文的重点，是讨论独立数独的个数怎么求。所谓独立数独的个数，就是在对称性操作下，依然不同的数独的个数。这里的对称性操作则是指的对于合法数独排布操作之后依旧是合法数独的变换。n=2 时的对称性主要有以下这些：
- 12 行交换，34 行交换，12 大行和 34 大行交换 （列同理）
- 横纵及两条对角线，4个镜面操作
- 90，180，270 旋转
- 填充数字 1，2，3，4 的 24 个置换操作
很明显这些操作构成一个非阿贝尔群。其中操作 4 和前三个操作是完全对易的。因此我们可以这样考虑，我们把第一子块固定为 1234 的数独解成为标准解。我们作用任意多次 1，2，3中的操作到一个标准解上，最后通过 4 中对应的置换操作，把第一子块换回 1234，那么此时又对应一个标准解。任何可以通过以上过程连接的标准解都对应了一个独立解。更严格的说，我们想找到独立数独的个数，在群论里就是要找到轨道数。所谓轨道就是群 G 中的元素作用在集合 X 上，那么 X 中元素通过群作用的相互变换会形成不交叉不空白的若干轨道，每个轨道里包含的 X 中的元素都可以通过群 G 中的元素作用互相得到。这里标准数独解构成集合 X，而 1，2，3，4条里提到的对称性构成了群 G。那么对称性下独立的标准解的个数（非标准解一定不独立，因为可从标准解置换 1234 得到），恰好就是对称群作用在标准解集合上的轨道数目。而这一轨道数目的求法，则是非常著名的 Burnside 引理4。该引理指出，群作用在某集合上的轨道的数目恰好等于每个群元素作用在该集合上不动点的数目的平均值。讲解该引理最浅显的例子，就是计算考虑对称性下，方形或立方体可能的染色情况的总数，该例子可参考上述 wiki 的链接，这里不再赘述。我们考虑的群 G，是保标准解集合 X 封闭的全部操作。
通过这种思路，我们需要求的东西包括群 G 的共轭类的数目及代表元素（每个共轭类的元素对应的 X 中的不动点数目相同），然后对每个类的代表元素求其不动点的数目。最后类元素数加权平均的不动点数即为考虑对称性后独立的数独的个数。下面展示的是实现这个思想的非常粗糙的代码，由于 n = 2 的数独问题规模很小，因此没做什么优化，只是示意一下（连矩阵 hash 都没做）。对于 n = 3 的问题，可以用专业的有限群计算软件5来计算群的阶数类数不动点等等。
def removearray(L, arr):
ind = 0
size = len(L)
while ind != size and not np.allclose(L[ind], arr):
ind += 1
if ind != size:
L.pop(ind)
def back_to_canonical(st0, first_grid=None):
st = st0.copy()
trs = {}
if not first_grid:
first_grid = [1, 2, 3, 4]
trs[st[0, 0]] = first_grid[0]
trs[st[0, 1]] = first_grid[1]
trs[st[1, 0]] = first_grid[2]
trs[st[1, 1]] = first_grid[3]
for i in range(4):
for j in range(4):
st[i, j] = trs[st0[i, j]]
return st
def r(i, st0):
st = st0.copy()
if i == 0:
return st
elif i == 1: # row exchange
st[0] = st0[1]
st[1] = st0[0]
elif i == 2: # row exchange
st[2] = st0[3]
st[3] = st0[2]
elif i == 3: # row block exchange
st[0] = st0[2]
st[1] = st0[3]
st[2] = st0[0]
st[3] = st0[1]
elif i == 4: # col exchange
st[:, 1] = st0[:, 0]
st[:, 0] = st0[:, 1]
elif i == 5: # col exchange
st[:, 2] = st0[:, 3]
st[:, 3] = st0[:, 2]
elif i == 6: # col block exchange
st[:, 0] = st0[:, 2]
st[:, 1] = st0[:, 3]
st[:, 2] = st0[:, 0]
st[:, 3] = st0[:, 1]
elif i == 7: # mirror on x
st[0] = st0[3]
st[1] = st0[2]
st[2] = st0[1]
st[3] = st0[0]
elif i == 8: # mirror on y
st[:, 0] = st0[:, 3]
st[:, 1] = st0[:, 2]
st[:, 2] = st0[:, 1]
st[:, 3] = st0[:, 0]
elif i == 9: # mirror on diagonal
st = np.transpose(st0)
elif i == 10: # mirror on another diagonal
st = np.rot90(st0, 2).T
elif i == 11: # rotation
st = np.rot90(st0)
elif i == 12: # rotation
st = np.rot90(st0, 2)
elif i == 13: # rotation
st = np.rot90(st0, 3)
return st
def apply_element(op, st0):
st = st0.copy()
for i in range(4):
for j in range(4):
st[i, j] = st0[int(op[i, j] // 4), int(op[i, j] % 4)]
return st
def apply_element_inverse(op, st0):
st = st0.copy()
for i in range(4):
for j in range(4):
st[int(op[i, j] // 4), int(op[i, j] % 4)] = st0[i, j]
return st
def find_all_elements(): # find all elements of symmetry group from generators defined by r()
st = np.arange(16).reshape([4, 4])
sts = [st]
q = [st]
while True:
for i in range(1, 14):
nst = r(i, q[0])
for oldst in sts:
if np.allclose(nst, oldst):
break
else:
sts.append(nst)
q.append(nst)
if len(q) > 1:
q = q[1:]
else:
break
return sts
def find_all_classes(elements): # find all conjugacy classes of the group
bc = np.arange(16).reshape([4, 4])
dq = {}
unclassfied = elements
while True:
e = unclassfied[0]
newno = len(dq) + 1
dq[newno] = [e]
for _, gs in dq.items():
for g in gs:
esame = apply_element_inverse(g, apply_element(e, apply_element(g, bc)))
for olde in dq[newno]:
if np.allclose(olde, esame):
break
else:
dq[newno].append(esame)
removearray(unclassfied, esame)
for g in unclassfied[1:]:
esame = apply_element_inverse(g, apply_element(e, apply_element(g, bc)))
for olde in dq[newno]:
if np.allclose(olde, esame):
break
else:
dq[newno].append(esame)
removearray(unclassfied, esame)
if len(unclassfied) > 1:
unclassfied = unclassfied[1:]
else:
break
return dq
def solutions():
r = []
st = np.zeros([4, 4])
st[0, 0] = 1
st[0, 1] = 2
st[1, 0] = 3
st[1, 1] = 4
return find_solution(0, 0, st, r)
sols = solutions()
def fixed_point(op):
r = 0
for sol in sols:
if np.allclose(back_to_canonical(apply_element(op, sol)), sol):
r += 1
return r
eles = find_all_elements()
clas = find_all_classes(eles)
print(len(eles), len(clas)) # 128 20
print([(k, len(v)) for k, v in clas.items()])
print(sum([len(v) * fixed_point(v[0]) for k, v in clas.items()]) / len(eles)) # 2
由此程序，可直接计算得到 n = 2 独立数独解的个数为 2。而且可以知道 n = 2 数独对称群的阶数是 128， 类数是 20。这两个独立解分别是 1234，3412，2143，4321 和 1234，3412，2341，4123。可以用以下代码验证这两个解确实无法被对称群相互变换。
def check_independence(isols, eles):
for s1 in isols:
for s2 in isols:
if not np.allclose(s1, s2):
for e in eles:
assert not np.allclose(back_to_canonical(apply_element(e, s1)), s2)
isols = [
np.array([[1, 2, 3, 4], [3, 4, 1, 2], [2, 1, 4, 3], [4, 3, 2, 1]]),
np.array([[1, 2, 3, 4], [3, 4, 1, 2], [2, 3, 4, 1], [4, 1, 2, 3]]),
]
check_independence(isols, eles)
事实上，对于 n=2 的情况，我们完全可以以 12 个标准解为图的顶点，以对称操作的连接为图的边，都穷举出来之后，在 BFS 搜索，查看 disconnected 图的数目，即可确定有几个独立解。但是本文的方法论都是对 n=3 的预热，可以看到利用 Burnside 引理的上述方法，可以解决 n=3 数独独立解数目的问题，但直接穷举是不行的。
关于用上述群论方法对 n = 3 数独独立解的计数，一些结果如下6。n=3 数独对称群，阶数为 3359232，共轭类数为 275，计算得到的独立解的个数为 5472730538。也就是说考虑了这么多的交换，置换，旋转，镜面对称性，独立解的个数还是几十亿的数量级，这也是比较违反我的直觉的。数独独立解的个数对应整数序列 A109741。
References
-
http://www.afjarvis.staff.shef.ac.uk/sudoku/，事实上，全部 n=3 数独的个数在 2003 年即被某论坛用户 QSCGZ 算出。参见 https://groups.google.com/forum/#!topic/rec.puzzles/A7pi7S12oFI， 该用户也同时给出了考虑对称性后的不同数独的数目估计，但似乎和真实值偏差较大 （似乎是真实值的一半，难道现在的精确结果漏算某个对称性？）。 ↩