10000个状态，全部找到最优解——我的魔表求解器和上帝数对上了
前文得到了魔表的 30 操作通解公式，包含 16 个自由变量。本文将其转化为 ℓ0 稀疏优化问题，用 Google OR-Tools CP-SAT 约束求解器求全局最优解。对一万个随机状态求解，100% 证得全局最优，步数分布与 cube20.org 穷举结果在 0.5% 精度内完全吻合。两种操作模型殊途同归，上帝数为 12。
魔表07：最少步数——把最优解交给约束求解器
回顾
在上一篇文章中，我们将魔表 \(14 \times 30\) 矩阵 \(\mathbf{M}\) 在 \(\mathbb{Z}_{12}\) 下做了高斯消元，得到了预变换法的核心公式：
其中 \(\mathbf{E}\)（\(14 \times 14\)）是行变换矩阵，\(\mathbf{P}\)（\(30 \times 30\)）是列置换矩阵，\(\mathbf{F}\)（\(14 \times 16\)）是自由列系数矩阵。\(\mathbf{t}\) 是 \(16\) 维自由变量，可取任意 \(0 \sim 11\) 的值。
如上一篇结尾所说，本篇的目标是：
在 \(12^{16}\) 种 \(\mathbf{t}\) 中，找到使 \(\|\mathbf{k}\|_0\) 最小的那个——也就是魔表的最优解。
一、问题重述
1.1 我们要最小化什么
令 \(\mathbf{k}_{\text{pivot}} = \mathbf{E}\mathbf{b} - \mathbf{F}\mathbf{t}\)（14 维），\(\mathbf{k}_{\text{free}} = \mathbf{t}\)（16 维）。经过列置换 \(\mathbf{P}\) 还原后，最终操作向量 \(\mathbf{k}\) 的非零分量数等于两者非零分量数之和。
记 \(\|\cdot\|_0\) 为 \(\ell_0\)范数（向量中非零分量的个数），我们的目标是：
换言之：在 \(12^{16}\) 种 \(\mathbf{t}\)（约 \(1.8 \times 10^{17}\) 种）中，找到让总操作数最少的那一个。
1.2 博客例子
继续用第二篇文章中的打乱状态。14 维简化状态为：
取负得 \(\mathbf{b} = -\mathbf{x}_{14} \bmod 12\)。通过预变换得 \(\mathbf{b}' = \mathbf{E}\mathbf{b}\)：
令 \(\mathbf{t} = \mathbf{0}\) 时得到特解：13 步（详见博客 06 第六节）。但这不是最优解——13 步中很可能有些操作是冗余的，可以通过选一个合适的 \(\mathbf{t} \neq \mathbf{0}\) 来消除。
二、F 矩阵的一个细节
在进入优化之前，有一个关于 \(\mathbf{F}\) 矩阵的性质值得注意。
2.1 实数域 vs 模 12
\(\mathbf{F}\) 是 \(14 \times 16\) 矩阵。在实数域 \(\mathbb{R}\) 上，\(\mathbf{F}\) 的秩是 \(14\)（满行秩）。但在 \(\mathbb{Z}_{12}\) 下，秩只有 \(13\)——比满行秩少 \(1\)。
这意味着 \(\mathbf{F}\) 在 \(\mathbb{Z}_{12}\) 下有一个非零左零向量：
满足 \(\mathbf{w}^T\mathbf{F} \equiv \mathbf{0} \pmod{12}\)（对所有 16 列）。
这也就是说，\(\mathbf{F}\) 矩阵的行2, 行6, 行12 线性相关：
2.2 这意味着什么
对任意 \(\mathbf{t}\)，左乘 \(\mathbf{w}^T\)：
\(\mathbf{w}^T\mathbf{E}\mathbf{b}\) 的值不随 \(\mathbf{t}\) 改变。对于博客例子：
由于结果 \(\mathbf{w}^T\mathbf{E}\mathbf{b} \neq 0\)，所以本例对于任意的 \(\mathbf{t}\) 来说，\(\mathbf{E}\mathbf{b} - \mathbf{F}\mathbf{t}\) 中位置 \(\{2, 6, 12\}\) 至少有一个不能变为 \(0\)——也就是说至少 \(1\) 个 pivot 操作是强制的。
需要说明的是，这个约束在第三节的 CP-SAT 建模中不需要显式添加，此条件被隐含满足。换句话说，\(\mathbf{w}\) 提供的是关于解空间结构的理论洞察，但在代码层面，它并不对应任何显式的约束。
这个信息对优化有方向性帮助，但本质上，我们还是需要系统地搜索最优 \(\mathbf{t}\)。
三、用 CP-SAT 建模
3.1 为什么是约束求解器
\(\ell_0\) 最小化在整数域上是 NP 难问题。直接枚举 \(12^{16}\) 种 \(\mathbf{t}\) 不可能。但现代约束规划（Constraint Programming, CP）求解器专门处理这类组合优化问题。
我们选择 Google 的 OR-Tools CP-SAT，一个成熟的有限域整数优化求解器。它内部实现了分支定界、冲突分析、域传播等算法，我们只需要把问题用数学约束表达出来，剩下的交给它。
3.2 变量与约束
引入三类变量：
| 变量 | 类型 | 含义 |
|---|---|---|
| \(x_0,\dots,x_{13}\) | 整数，\(0\sim11\) | \(\mathbf{k}_{\text{pivot}}\) 的 14 个分量 |
| \(t_0,\dots,t_{15}\) | 整数，\(0\sim11\) | \(\mathbf{t}\) 的 16 个分量 |
| \(p_0,\dots,p_{13},\;q_0,\dots,q_{15}\) | 布尔，\((0/1)\) | 指示对应的 \(x\) 或 \(t\) 是否为非零 |
约束 1：模等式的显式形式。 对于每行 \(i\)，引入辅助整数变量 \(y_i\)（代表“借 12”的次数），将模等式写为普通线性等式：
\(y_i\) 的值域并非任意。从 \(\sum_j F_{ij} \cdot t_j\) 的最大值可以推算 \(y_i\) 的精确上界（约 \(0 \sim 100\) 之间，取决于 \(F\) 各行的稀疏度）。等式两边都非负，避免了模运算引入的符号复杂性。
为什么不用
AddModuloEquality？ OR-Tools 提供了内置的模等式约束。但在实际测试中，我们发现它对大表达式（涉及 16 个 \(t\) 变量）的内部分解存在域传播问题——会将某些数学上正确的解错误地判定为不可行。显式等式x + Ft = Eb + 12y绕开了这个问题，并且数学上完全等价。
约束 2：指示变量的含义。 用半具体化约束（OnlyEnforceIf）建立 \(p_i \leftrightarrow (x_i \neq 0)\) 和 \(q_j \leftrightarrow (t_j \neq 0)\) 的对应关系：
- \(p_i = 0 \;\Longrightarrow\; x_i = 0\)
- \(p_i = 1 \;\Longrightarrow\; x_i \neq 0\)
对 \(q_j\) 和 \(t_j\) 同理。
目标函数：
整个模型 74 个变量、74 条约束（直接读取 CP-SAT 内部 proto 的统计值，含 14 个整数 \(x\)、16 个整数 \(t\)、14 个整数 \(y\)、14 个布尔 \(p\)、16 个布尔 \(q\) 以及等价的线性等式、指示约束和目标函数）。CP-SAT 可以高效处理。
3.3 关键代码
def solve_optimal(b_prime, F):
"""返回最小步数和对应的 x, t"""
model = cp_model.CpModel()
# 变量
t = [model.NewIntVar(0, 11, f't{j}') for j in range(16)]
x = [model.NewIntVar(0, 11, f'x{i}') for i in range(14)]
y = [model.NewIntVar(0, y_max[i], f'y{i}') for i in range(14)]
p = [model.NewBoolVar(f'p{i}') for i in range(14)]
q = [model.NewBoolVar(f'q{j}') for j in range(16)]
# 约束 1：x + Ft = b' + 12y
for i in range(14):
ft = sum(F[i][j] * t[j] for j in range(16))
model.Add(x[i] + ft == b_prime[i] + 12 * y[i])
# 约束 2：p ↔ x≠0, q ↔ t≠0
for i in range(14):
model.Add(x[i] == 0).OnlyEnforceIf(p[i].Not())
model.Add(x[i] != 0).OnlyEnforceIf(p[i])
for j in range(16):
model.Add(t[j] == 0).OnlyEnforceIf(q[j].Not())
model.Add(t[j] != 0).OnlyEnforceIf(q[j])
# 目标
model.Minimize(sum(p) + sum(q))
# 求解
solver = cp_model.CpSolver()
status = solver.Solve(model)
x_opt = [solver.Value(x[i]) for i in range(14)]
t_opt = [solver.Value(t[j]) for j in range(16)]
steps = int(solver.ObjectiveValue())
return steps, x_opt, t_opt
点击查看求解器完整 Python 代码
"""
魔表最少步数求解器
==================
更新时间: 2026-07-08 21:25
作者: H_Elden (https://www.cnblogs.com/h-elden)
==================
基于 5 个原型 + P/Q 联动，直接生成 30 种基本操作。
核心改进（相比 6 原型旧方案）：
1. 原型从 6 个减为 5 个：v3 = R(v2)，由 P/Q 联动推导，无需单独测量。
2. 不再"先生成 64 列再去重"，而是利用 P/Q 联动直接挑选 30 个代表操作。
3. 原型向量直接用 14 维（去掉背面 4 角的冗余指针），无需后续删行。
操作记法：OptClock 记法
每种操作用 5 个字符表示，如 UUDD u3'
- 4 大写字母 = 4 个 pin 状态，顺序：左上(UL) 右上(UR) 左下(DL) 右下(DR)
U = 弹起 (向你), D = 按下 (远离你)
- 1 小写字母 = 转动哪个拨轮
u = 转 U pin 旁的拨轮 (F 轮), d = 转 D pin 旁的拨轮 (B 轮)
- 转动量：字母本身 = 顺时针 1 格；数字 = 顺时针 n 格；撇号 = 逆时针
u3 = 顺时针 3 格
u3' = 逆时针 3 格 (= 顺时针 9 格)
u' = 逆时针 1 格 (= 顺时针 11 格)
14 维坐标说明：
18 维 = [正面 3×3 (9), 反面 3×3 (9)]
背面 4 角 (B1,B3,B7,B9) 由不变量自动确定：B1=-F3, B3=-F1, B7=-F9, B9=-F7
14 维 = [F1..F9, B2, B4, B5, B6, B8] ← 去掉背面 4 角
数学模型：M_{14×30} · k ≡ -x (mod 12)
- M：14×30 操作矩阵
- k：30 维操作向量（每种操作的转动格数，0~11）
- x：14 维指针状态（需要复原的打乱状态）
求解流程：
5 原型 → C/D 生成 16 种 UL 向量 → P/Q 联动直接选 30 操作
→ 构建 M → Z_12 高斯消元得 E, F → CP-SAT 求 L0 最稀疏解
"""
import time
import numpy as np
from ortools.sat.python import cp_model
# ============================================================
# 第一部分：14↔18 维转换与对称变换算子 R, C, D
# ============================================================
# 14 维 = [F1..F9, B2, B4, B5, B6, B8]
# 背面 4 角由不变量补全：B1=-F3, B3=-F1, B7=-F9, B9=-F7
# 18 维中保留的行索引（去掉 B1=9, B3=11, B7=15, B9=17）
_ROW_14 = [0, 1, 2, 3, 4, 5, 6, 7, 8, 10, 12, 13, 14, 16]
def to_18(v14):
"""14 维 → 18 维（用不变量补全背面 4 角）"""
v18 = np.zeros(18, dtype=int)
v18[_ROW_14] = v14 % 12
v18[9] = (-v18[2]) % 12 # B1 = -F3
v18[11] = (-v18[0]) % 12 # B3 = -F1
v18[15] = (-v18[8]) % 12 # B7 = -F9
v18[17] = (-v18[6]) % 12 # B9 = -F7
return v18
def to_14(v18):
"""18 维 → 14 维（去掉背面 4 角）"""
return v18[_ROW_14]
def to_matrices(v18):
"""18 维向量 → (正面 3×3, 反面 3×3)"""
return v18[:9].reshape(3, 3), v18[9:].reshape(3, 3)
def to_vector(F, B):
"""(正面 3×3, 反面 3×3) → 18 维向量"""
return np.concatenate([F.flatten(), B.flatten()])
def R_op(v14, k=1):
"""旋转 R^k：正面顺时针 k×90°，反面逆时针 k×90°（14 维输入输出）"""
v18 = to_18(v14)
F, B = to_matrices(v18)
return to_14(to_vector(np.rot90(F, k=-int(k)), np.rot90(B, k=int(k))))
def C_op(v14):
"""正反互补：C(F,B) = H(-B, -F)，H 为水平镜像（14 维输入输出）"""
v18 = to_18(v14)
F, B = to_matrices(v18)
return to_14(to_vector(np.fliplr(-B), np.fliplr(-F)))
def D_op(v14):
"""对角镜像：正面转置，反面反对角转置（14 维输入输出）"""
v18 = to_18(v14)
F, B = to_matrices(v18)
return to_14(to_vector(F.T, np.rot90(B, k=2).T))
# ============================================================
# 第二部分：5 个原型效果向量（14 维）
# ============================================================
# 固定拨轮为 UL、顺时针转 1 格，手动测量的 5 个基本原型。
# 第 6 个原型 v3 已被证明等于 R(v2)（P/Q 联动），不再单独测量。
#
# OptClock 记法：4 大写字母 = pin 状态 (UL UR DL DR)
# U = 弹起, D = 按下; 小写 u = 转 U pin 旁拨轮, d = 转 D pin 旁拨轮
# 所有原型拨轮固定为 UL，故小写字母取决于 UL 的 pin 状态。
v0 = np.array([1,1,1,1,1,1,1,1,1, 0,0,0,0,0], dtype=int) # UUUU u 无按钮按下
v1 = np.array([1,0,0,0,0,0,0,0,0, 11,0,11,11,0], dtype=int) # DUUU d 只按 UL
v2 = np.array([1,1,0,1,1,1,1,1,1, 0,0,0,0,0], dtype=int) # UDUU u 只按 UR
v4 = np.array([1,0,1,0,0,0,0,0,0, 11,11,11,11,0], dtype=int) # DDUU d 按 UL+UR
v5 = np.array([1,0,0,0,0,0,0,0,1, 11,11,11,11,11], dtype=int) # DUUD d 按 UL+DR
# v3 = R(v2)：R 将 UL 拨轮的效果旋转到 UR 拨轮，
# 而 UR 与 UL 同在 Q 组（P/Q 联动），效果等价。
# 对应操作 UUUD u：只按 DR，转 U pin 旁拨轮。
v3 = R_op(v2) % 12 # ← 推导得出，无需测量
# ============================================================
# 第三部分：用 C, D 从 5 原型生成 16 种按钮状态（UL 拨轮）
# ============================================================
# 4 个按钮共 2^4 = 16 种按下/弹起组合。
# 每个原型与其补集（C 变换）交替排列，D 变换补上对角线镜像。
# label 使用 OptClock 记法：pin 串 + 小写字母。
# 小写 u/d 由 UL 的 pin 状态决定（U→u, D→d）。
ul_vecs, ul_lbls = [], []
def add(label, vec):
ul_vecs.append(np.array(vec, dtype=int) % 12)
ul_lbls.append(label)
add('UUUU u', v0) # 无按钮
add('DDDD d', C_op(v0)) # 全按下
add('DUUU d', v1) # 只按 UL
add('UDDD u', C_op(v1)) # 按 UR+DL+DR
add('UDUU u', v2) # 只按 UR
add('DUDD d', C_op(v2)) # 按 UL+DL+DR
add('UUDU u', D_op(v2)) # 只按 DL
add('DDUD d', C_op(D_op(v2))) # 按 UL+UR+DR
add('UUUD u', v3) # 只按 DR (= R(v2))
add('DDDU d', C_op(v3)) # 按 UL+UR+DL
add('DDUU d', v4) # 按 UL+UR
add('UUDD u', C_op(v4)) # 按 DL+DR
add('DUDU d', D_op(v4)) # 按 UL+DL
add('UDUD u', C_op(D_op(v4))) # 按 UR+DR
add('DUUD d', v5) # 按 UL+DR
add('UDDU u', C_op(v5)) # 按 UR+DL
# ============================================================
# 第四部分：P/Q 联动直接生成 30 种基本操作
# ============================================================
# P/Q 联动：在固定按钮状态下，四个拨轮分成两个等价组——
# B 轮（按下按钮对应的拨轮，P 组）：组内任意拨轮效果相同 → OptClock 小写 d
# F 轮（弹起按钮对应的拨轮，Q 组）：组内任意拨轮效果相同 → OptClock 小写 u
#
# 因此每种按钮状态只需取至多 2 个代表拨轮（1 个 B 轮 + 1 个 F 轮）：
# UUUU u（全弹起）→ 只有 F 轮，4 个全等 → 1 种操作
# DDDD d（全按下）→ 只有 B 轮，4 个全等 → 1 种操作
# 其余 14 种 → 1 个 F 轮 + 1 个 B 轮 → 2 种操作
# 合计：1 + 1 + 14×2 = 30 种
#
# 不再生成 64 列再去重，直接挑选代表。
# OptClock pin 串 → 4-bit 按钮编码
# pin 串顺序：UL UR DL DR；bit 编码：bit0=UL, bit1=UR, bit2=DR, bit3=DL
def pin_to_bits(pin):
b = 0
if pin[0] == 'D': b |= 1 # UL
if pin[1] == 'D': b |= 2 # UR
if pin[3] == 'D': b |= 4 # DR
if pin[2] == 'D': b |= 8 # DL
return b
def rotate_bits(bits, k):
"""4-bit 左旋转 k 位 (k>0 = R 方向: UL→UR→DR→DL)"""
k = k % 4
return ((bits << k) | (bits >> (4 - k))) & 0xF
wheels = ['UL', 'UR', 'DR', 'DL'] # R 旋转方向
# UL 向量查找表：pin bits → 14 维向量
lookup = {pin_to_bits(lbl[:4]): v for lbl, v in zip(ul_lbls, ul_vecs)}
# 旋转扩展公式：e_{S, R^k(UL)} = R^k(e_{R^{-k}(S), UL})
def effect_vector(S_bits, wheel_idx):
"""计算按钮集 S 在拨轮 wheel_idx 下的效果向量 (14 维, mod 12)"""
k = wheel_idx
S0_bits = rotate_bits(S_bits, -k) # R^{-k}(S)
vec14 = lookup[S0_bits] # E(R^{-k}(S), UL) 14 维
return R_op(vec14, k) % 12 # R^k(...) 14 维
# 判断拨轮在当前按钮状态下是 B 轮还是 F 轮
def is_b_wheel(pin, wheel):
"""拨轮 wheel 对应的 pin 是否按下 (D)"""
idx = {'UL': 0, 'UR': 1, 'DL': 2, 'DR': 3}
return pin[idx[wheel]] == 'D'
# 直接生成 30 种操作：每种按钮状态取各组首个代表拨轮
# unique_info: [(14维向量, OptClock label), ...]
# label = pin串 + 空格 + 小写字母 (u 或 d)
unique_info = []
for label, _ in zip(ul_lbls, ul_vecs):
pin = label[:4]
S_bits = pin_to_bits(pin)
seen_B = False # 是否已取 B 轮代表
seen_F = False # 是否已取 F 轮代表
for k, w in enumerate(wheels):
if is_b_wheel(pin, w):
if not seen_B:
seen_B = True
unique_info.append((effect_vector(S_bits, k), f"{pin} d"))
else:
if not seen_F:
seen_F = True
unique_info.append((effect_vector(S_bits, k), f"{pin} u"))
assert len(unique_info) == 30, f"预期 30 种操作，实际 {len(unique_info)}"
# 30 种操作的 OptClock label 列表
opt_labels = [info[1] for info in unique_info]
# ============================================================
# 第五部分：构建 14×30 矩阵 M
# ============================================================
# 原型向量已经是 14 维，M 矩阵直接构建，无需删行。
M = np.zeros((14, 30), dtype=int)
for c, (vec, _) in enumerate(unique_info):
M[:, c] = vec
# ============================================================
# 第六部分：Z_12 高斯消元 → E, F, col_perm
# ============================================================
# 将 M 消元为 [I | F] 形式，同时记录行变换 E 和列置换 col_perm。
# 模 12 下 {1,5,7,11} 自逆，可直接作为主元归一用的乘法逆元。
m, n = M.shape # 14, 30
E_mat = np.eye(m, dtype=int)
A_work = M.copy() % 12
col_perm = list(range(n)) # col_perm[新位置] = 原始列号
row = 0
for col in range(n):
if row >= m:
break
# 在剩余子矩阵中找非零主元
found_r, found_c = None, None
for c2 in range(col, n):
for r in range(row, m):
if A_work[r, c2] != 0:
found_r, found_c = r, c2
break
if found_r is not None:
break
if found_r is None:
continue
# 列交换
if found_c != col:
A_work[:, [col, found_c]] = A_work[:, [found_c, col]]
col_perm[col], col_perm[found_c] = col_perm[found_c], col_perm[col]
# 行交换
if found_r != row:
A_work[[row, found_r]] = A_work[[found_r, row]]
E_mat[[row, found_r]] = E_mat[[found_r, row]]
# 主元归一
pivot = A_work[row, col]
inv = pow(int(pivot), -1, 12)
A_work[row] = (A_work[row] * inv) % 12
E_mat[row] = (E_mat[row] * inv) % 12
# 消去其他行
for r in range(m):
if r != row and A_work[r, col] != 0:
factor = A_work[r, col]
A_work[r] = (A_work[r] - factor * A_work[row]) % 12
E_mat[r] = (E_mat[r] - factor * E_mat[row]) % 12
row += 1
# 自由列系数矩阵 F (14×16)
F_mat = A_work[:14, 14:] % 12
F_list = [[int(F_mat[i, j]) for j in range(16)] for i in range(14)]
# y 变量的紧界（减少搜索空间）
_max_ft = [sum(F_list[i][j] * 11 for j in range(16)) for i in range(14)]
_y_max = [int((11 + mf) / 12) for mf in _max_ft]
# ============================================================
# 第七部分：CP-SAT 最少步数求解函数
# ============================================================
# 将通解 k = P⁻¹·[x; t] 代入 M·k ≡ b (mod 12) 后，
# 最少步数 = min ||k||_0 = min (非零 x 个数 + 非零 t 个数)
# 用布尔指示变量 p_i, q_j 表示 x_i, t_j 是否非零，最小化 sum(p) + sum(q)。
def solve_optimal(Eb_list, F_list, y_max):
"""CP-SAT 求解最少步数
输入：Eb_list = E·b (mod 12) 的 14 个分量
F_list = 自由列系数矩阵 (14×16)
y_max = y 变量的上界
返回：(steps, x_opt, t_opt, status_name) 或 (None, None, None, status_name)"""
model = cp_model.CpModel()
# 变量
t = [model.NewIntVar(0, 11, f't{j}') for j in range(16)]
x = [model.NewIntVar(0, 11, f'x{i}') for i in range(14)]
y = [model.NewIntVar(0, y_max[i], f'y{i}') for i in range(14)]
p = [model.NewBoolVar(f'p{i}') for i in range(14)]
q = [model.NewBoolVar(f'q{j}') for j in range(16)]
# 约束 1：x + Ft = b' + 12y
for i in range(14):
ft = sum(F_list[i][j] * t[j] for j in range(16))
model.Add(x[i] + ft == Eb_list[i] + 12 * y[i])
# 约束 2：p ↔ x≠0, q ↔ t≠0
for i in range(14):
model.Add(x[i] == 0).OnlyEnforceIf(p[i].Not())
model.Add(x[i] != 0).OnlyEnforceIf(p[i])
for j in range(16):
model.Add(t[j] == 0).OnlyEnforceIf(q[j].Not())
model.Add(t[j] != 0).OnlyEnforceIf(q[j])
# 目标：最少步数 = 非零操作数
model.Minimize(sum(p) + sum(q))
# 求解
solver = cp_model.CpSolver()
status = solver.Solve(model)
if status in [cp_model.OPTIMAL, cp_model.FEASIBLE]:
x_opt = [solver.Value(x[i]) for i in range(14)]
t_opt = [solver.Value(t[j]) for j in range(16)]
steps = int(solver.ObjectiveValue())
return steps, x_opt, t_opt, solver.StatusName(status)
return None, None, None, solver.StatusName(status)
def solve_state(x14):
"""求解初始状态 x14 的最少步数复原方案
输入：x14 = 14 维指针状态向量 (0~11)
返回：(steps, k_full, status) 或 (None, None, status)
- steps: 最少步数
- k_full: 30 维操作向量 (每种操作的转动格数 0~11)
- status: 求解器状态字符串"""
b14 = (-x14) % 12
Eb = (E_mat @ b14) % 12
Eb_list = [int(Eb[i]) for i in range(14)]
steps, x_opt, t_opt, status = solve_optimal(Eb_list, F_list, _y_max)
if steps is None:
return None, None, status
# 还原 30 维解向量 k
k_perm = np.zeros(30, dtype=int)
for i in range(14):
k_perm[i] = x_opt[i]
for j in range(16):
k_perm[14 + j] = t_opt[j]
k_full = np.zeros(30, dtype=int)
for np_idx in range(30):
orig = col_perm[np_idx]
k_full[orig] = k_perm[np_idx]
return steps, k_full, status
# ============================================================
# 第八部分：交互式求解
# ============================================================
def optclock_move(label, k):
"""OptClock label + 系数 k (0~11) → 完整转动记法；k=0 返回 None"""
k = k % 12
if k == 0:
return None
if k <= 6:
return label if k == 1 else f"{label}{k}"
else:
ccw = 12 - k
return f"{label}'" if ccw == 1 else f"{label}{ccw}'"
row_names_14 = ['F1','F2','F3','F4','F5','F6','F7','F8','F9','B2','B4','B5','B6','B8']
def main():
"""交互式求解入口"""
print("=" * 60)
print(" 魔表最少步数求解器 (5 原型 · 30 操作 · 14 维)")
print("=" * 60)
print()
print("OptClock 记法：UUDD u3' = pin状态(UUDD) + 转U轮逆时针3格")
print(" U=弹起 D=按下 | u=转U轮 d=转D轮 | 撇号=逆时针")
print()
# 输出全部 30 种基本操作
# print("30 种基本操作：")
# print(f" {'#':>2s} {'操作':<8s}")
# print(" " + "-" * 20)
# for i, label in enumerate(opt_labels):
# print(f" {i+1:>2d} {label:<8s}")
# print()
print("输入格式：14 个 0~11 的数字，空格分隔")
print(f"顺序：{' '.join(row_names_14)}")
print("示例：8 2 9 4 1 0 2 4 10 9 3 11 6 3")
print("输入 'q' 退出")
print()
while True:
user_input = input("请输入 14 个数字 > ").strip()
if user_input.lower() == 'q':
break
try:
nums = [int(x) for x in user_input.split()]
except ValueError:
print("错误：请输入数字，空格分隔")
continue
if len(nums) != 14:
print(f"错误：需要 14 个数字，你输入了 {len(nums)} 个")
continue
if any(n < 0 or n > 11 for n in nums):
print("错误：每个数字必须在 0~11 之间")
continue
x14 = np.array(nums, dtype=int)
# 求解
print("求解中...", end=" ", flush=True)
t0 = time.perf_counter()
steps, k_full, status = solve_state(x14)
elapsed = time.perf_counter() - t0
print(f"状态: {status}，耗时 {elapsed:.2f}s")
if steps is not None:
print(f"最少步数: {steps}")
# 验证 M·k ≡ b (mod 12)
b14 = (-x14) % 12
result = (M @ k_full) % 12
ok = np.all(result == b14)
print(f"验证: {'通过' if ok else '失败'}")
print()
print("操作方案:")
moves = []
for op_idx in range(30):
k = int(k_full[op_idx])
if k != 0:
move = optclock_move(opt_labels[op_idx], k)
moves.append(move)
print(f" {len(moves):2d}. {move}")
print()
print("紧凑格式:")
print(' '.join(moves))
else:
print("求解失败")
print()
if __name__ == '__main__':
main()
四、博客例子：从 13 步到 9 步
将博客例子的 \(\mathbf{b}'\) 输入 CP-SAT 模型：
- t = 0 特解：13 步
- CP-SAT 最优解：9 步，耗时约 0.14 秒，状态 OPTIMAL（已证明全局最优）
4.1 解的构成
CP-SAT 找到的 \(\mathbf{x}\) 和 \(\mathbf{t}\) 分别为：
\(\mathbf{x}\) 中有 \(5\) 个非零元，\(\mathbf{t}\) 中有 \(4\) 个非零元，合计 9 步。比 \(\mathbf{t} = \mathbf{0}\) 时的 13 步减少了 \(4\) 步。
将 \(\mathbf{x}, \mathbf{t}\) 按置换 \(\mathbf{P}\) 还原为 30 维的 \(\mathbf{k}\)，然后从 30 种操作中提取非零分量：
| 步骤 | 按钮状态 | 拨轮 | \(k_i\) | 动作 |
|---|---|---|---|---|
| 1 | 只按 UL | UL | 8 | 逆时针 4 格 |
| 2 | 只按 UL | UR | 9 | 逆时针 3 格 |
| 3 | 按 UL, DL, DR（不按 UR） | UR | 3 | 顺时针 3 格 |
| 4 | 按 UL, UR, DR（不按 DL） | DL | 11 | 逆时针 1 格 |
| 5 | 按 UL, UR | DR | 2 | 顺时针 2 格 |
| 6 | 按 DL, DR | DR | 2 | 顺时针 2 格 |
| 7 | 按 UR, DR | UR | 5 | 顺时针 5 格 |
| 8 | 按 UL, DR | UL | 8 | 逆时针 4 格 |
| 9 | 按 UL, DR | UR | 10 | 逆时针 2 格 |
代入验证：\(\mathbf{M}\mathbf{k} \equiv \mathbf{b} \pmod{12}\) ✓。9 步可解此状态。
4.2 交叉验证
为确认这个 9 步解确实是全局最优，我们用一位独立开发者公开发布的最优求解器 OptClock 对同一状态进行求解。该求解器也返回了 9 步，转换为我们的 30 操作格式后，非零分量完全一致。其输出内容如下：
9 moves is optimal:
DUDD u3 DDUD u' DUUU d4' DDUU u2 UUDD d2 UDUD d5 DUUU u3' DUUD u2' DUUD d4'
Optimal solution found in 0.46 seconds.
两个完全独立的实现得到了相同的 9 步解，互相印证。
五、大规模验证：与 cube20.org 对比
5.1 为什么要做大规模验证
我们 30 种操作模型得到的"最优步数"和 cube20.org 基于 15 个群生成元得到的“上帝数”，是否指向同一个分布？
这是检验我们模型正确性的终极测试。如果是，那么两种不同的操作模型就是等价的群表示。
5.2 实验设计
- 随机生成 \(10{,}000\) 个合法状态（14 个独立坐标，每个 \(0 \sim 11\) 均匀随机）
- 固定随机种子
np.random.seed(42)以便复现 - 每题限时 \(30\) 秒，用 CP-SAT 求全局最优
- 需要 OPTIMAL 状态（已证明最优），不接受 FEASIBLE（仅找到解但未证明最优）
5.3 结果
在 i9-12900H 平台上，10,000 题总计耗时约 \(31\) 分钟：
总样本数: 10000
已证明最优: 10000 (100.0%)
超时/未确定: 0 (0.0%)
平均耗时: 185ms/题
最慢耗时: 1705ms
无一超时。全部 10,000 个状态都证得了全局最优解。
5.4 步数分布
| 步数 | 样本数 | 样本占比 | cube20 总体占比 | 偏差 |
|---|---|---|---|---|
| 6 | 5 | 0.05% | 0.03% | +0.02% |
| 7 | 57 | 0.57% | 0.67% | −0.10% |
| 8 | 795 | 7.95% | 7.91% | +0.04% |
| 9 | 4146 | 41.46% | 41.13% | +0.33% |
| 10 | 4759 | 47.59% | 47.76% | −0.17% |
| 11 | 238 | 2.38% | 2.48% | −0.10% |
（距离 0~5 和 12 在总体中占比极低（合计 \(< 0.001\%\)），在 \(10{,}000\) 样本中不出现在预期之内。）
- 样本加权平均：9.43
- 总体加权平均（cube20）：9.43
由下图对比饼状图也可以直观看出，我们抽取 10000 个样本测试的结果，和总体分布数据偏差极小。
5.5 偏差有多小？
\(10{,}000\) 个随机样本中，步数 \(d\) 的样本比例 \(\hat{p}\) 近似服从正态分布，标准误：
以步数 9 为例（\(\hat{p} \approx 41.46\%\)，\(n = 10{,}000\)）：
\(95\%\) 置信区间为 \(\hat{p} \pm 2 \times SE \approx \pm 1.0\%\)。
上表中全部偏差都在 \(\pm 0.33\%\) 以内，远小于置信区间宽度。换言之：
在 \(0.5\%\) 的统计精度下，两种操作模型的最优步数分布完全一致。
这个结果还有另一层含义。cube20 的分布是用 BFS 穷举全部 \(12^{14}\) 个状态得到的"真实值"。我们的样本占比与真实分布高度吻合，说明我们 CP-SAT 求出的每一个最优步数都经得起统计检验。如果我们的算法存在系统性错误（比如经常返回次优解），分布就会整体偏大，根本不可能在 \(0.5\%\) 精度上对上。换言之，独立于任何求解器之外、通过统计一致性获得的证据，本身就验证了我们算法的正确性。
5.6 cube20 距离12 状态验证
除了随机状态，我们还从 cube20.org 公开发布的 dist12.txt（包含全部 \(39{,}248\) 个距离为 \(12\) 的“最远状态”）中提取了前 \(100\) 个，用我们的模型逐一求解。
距离 12 的“最远状态”比随机状态难得多——证明“不存在 11 步解”需要更深的搜索树。这 100 个状态平均每题 7.6 秒，全部 \(100\) 个总计 \(762\) 秒（约 \(13\) 分钟）。
全部 \(100\) 个状态都得到了 \(12\) 步解，异常数为 \(0\)。与 cube20 的上帝数结论完全吻合。
测试环境：Intel Core i9-12900H，16 GB RAM，Windows 11，Python 3.11.4，OR-Tools CP-SAT 9.15。CP-SAT 默认以单线程搜索为主，同一题目在不同 CPU 上的耗时主要取决于单核性能。
六、讨论：两种模型，一个分布
6.1 我们做了什么
用 \(30\) 种（button_set, wheel）操作构建 \(14 \times 30\) 矩阵 \(\mathbf{M}\)，在 \(\mathbb{Z}_{12}\) 下消元得到 \(\mathbf{E}, \mathbf{P}, \mathbf{F}\) 三个预计算矩阵，将最少步数求解转化为一个约束优化问题，交给 CP-SAT 求解器进行计算。
6.2 cube20 做了什么
用 \(15\) 种群生成元（touch pattern）在 \(9\) 维 coset 坐标上做双向 BFS 搜索，穷举全部 \(12^{14}\) 个状态，确定了上帝数 \(= 12\) 和完整的距离分布。
6.3 为什么分布一致
操作集不一样（\(30\) vs \(15\)），坐标表示不一样（\(14\) 维 vs \(9\) 维 coset），求解方法不一样（CP-SAT vs BFS），但最优步数分布完美吻合。
原因在于：尽管操作的表象不同，它们生成的群是相同的——\(\mathbb{Z}_{12}^{14}\)。两组操作互为等价生成元。上帝数是群的性质，不随生成元的选取而改变。
七、系列终章：回顾与总结
七篇博客，从零开始探索一个看似简单的小玩具（魔表 / Rubik's Clock），最终抵达了它的数学核心。
- 魔表01——建立模 12 运算基础，发现操作可交换、效果可线性叠加
- 魔表02——将全部 64 种操作用线性方程组建模，得到 \(18 \times 64\) 矩阵
- 魔表03——利用旋转 R、补 C、对角镜像 D 三种对称性，将 64 种操作归约为 6 个原型
- 魔表04——用代码从 6 个原型生成完整的 64 列矩阵 A
- 魔表05——解剖矩阵：发现 P/Q 联动、不变量、空间分解，将模型简化为 \(14 \times 30\)
- 魔表06——在 \(\mathbb{Z}_{12}\) 下直接消元，得到预变换法和包含 16 个自由变量的通解公式
- 本文——将 L0 稀疏优化建模为约束规划问题，用 CP-SAT 求全局最优，大规模验证与 cube20 结论一致
核心脉络：
数学建模是核心——矩阵 \(\mathbf{M}\)、行变换 \(\mathbf{E}\)、列置换 \(\mathbf{P}\)、通解公式，这些都来自对魔表结构的理解。CP-SAT 只是最后一公里的工具——它帮我们找到了最优解，但真正让这一切成为可能的是前面六篇博文打下的数学基础。
参考资料
Tom Rokicki，Rubik's Clock God's Number is 12，cube20.org.
OR-Tools CP-SAT Solver，Google.
链接：https://developers.google.com/optimization/cp/cp_solver
OptClock: an optimal Rubik's Clock solver，SpeedSolving 论坛.
链接：https://www.speedsolving.com/threads/optclock-optimal-rubiks-clock-solver.47747/
NP-hard 问题，百度百科.
浙公网安备 33010602011771号