蜂窝数独的SAT建模与C++求解：从CNF到DPLL

简介华中科技大学2022级程序设计综合课程设计任务一——基于SAT的蜂窝数独游戏求解程序面向计算机相关专业在校生及课程设计、毕设选题者。项目以DPLL算法为核心完整实现SAT求解器与蜂窝数独到CNF公式的转换涵盖CNF构建、答案处理、菜单交互等模块代码结构清晰、可直接运行。压缩包共11个文件以7个C源码文件为主体分别承担DPLL核心求解、CNF子句生成、谜题输出、耗时统计及交互菜单等功能辅以.h头文件、README说明、.gitignore及实验报告PDF整体体积仅1.23MB。实验报告包含问题分析、算法设计与测试结果适合作为课程设计参考或二次开发基础。该资源已有148人学习代码经测试通过答辩平均分96分对理解SAT求解、数独建模及C工程实践有较高参考价值也适合用于课程设计答辩演示与算法进阶练习。1. 蜂窝数独的SAT解法难的不是SAT而是网格建模把“蜂窝数独求解程序”做成课程设计最容易走偏的地方是一上来就写递归回溯填格子数字填不下去就撤销重来代码写了两三百行最后只对固定一张图有效。SAT 思路反着来先把“每格填什么数”表达成布尔变量再把“相邻不能相同”“同一轴线不能相同”写成 CNF 子句求解过程交给通用 SAT 算法。这个题目最值钱的部分其实是后者蜂窝网格的坐标、邻居关系、候选数字和布尔变量 id 一一对应规则怎么写决定了求解器是几百毫秒出结果还是跑到天亮。下面按一套能直接提交的 C 方案展开覆盖从 CNF 编码、六边形网格生成到 DPLL 实现的完整链路。2. 建模先行SAT变量定义与CNF子句生成2.1 变量编号一个格子一个候选值对应一个布尔变量蜂窝数独的 SAT 编码不是按“格”建变量而是按“格 × 候选值”建。设网格共有 N 个格子每个格子可填 1..K则变量总数是 N*K。变量编号直接取cellId*K (value-1)。比如第 5 格填 3候选数 K9就对应变量5*92。这样做的好处是约束翻译非常机械任意一条规则都能写成一组变量的析取或两两互斥而不用处理“一格多值”的特殊逻辑。变量编号的核心难点在于 cellId 必须稳定。不要用二维坐标直接做索引六边形网格用二维数组存会有大量空洞约束生成时遍历空洞会产生错误的邻居。我一般会先建一个vectorCell把有意义的格子压进去再用 compact 编号做映射。后续换网格形状或换候选数只改网格生成部分变量编号层完全不动。struct Cell { int q, r, s; // 轴向坐标s -q-r int id; // 0..N-1 }; int varId(int cellId, int value) { // value: 1..K return cellId * K (value - 1); // 0..N*K-1 } int litId(int cellId, int value) { // 1-based 文字0 保留给未赋值 return varId(cellId, value) 1; }varId返回 0-based 下标供数组索引用真正写进 CNF 子句的文字用litId正数代表该变量为真负数代表为假。这个区分必须从第一行代码就确立否则后面传播和回溯阶段会出现 0 号文字真假混淆的隐蔽 bug。2.2 每个格子必须且只能填一个值“必填一个值”对应一个长度为 K 的析取子句vectorint clause; for (int v 1; v K; v) clause.push_back(litId(cellId, v)); clauses.push_back(clause);“最多一个值”用两两互斥写最直观for (int v1 1; v1 K; v1) for (int v2 v1 1; v2 K; v2) { clauses.push_back({-litId(cellId, v1), -litId(cellId, v2)}); }候选数 K9 时每个格子产生 36 条二元子句19 格总共 684 条。对课程设计规模没问题如果题目把候选数字扩大到 16 或 25这里会明显膨胀到那时换成顺序编码sequential encoding加辅助变量把每格“至多一”的子句数从 O(K²) 降到 O(K)。实验报告里写清楚这个取舍本身是加分项。约束类型CNF 形式数量级每格至少一值(x1 ∨ x2 ∨ … ∨ xK)N每格至多一值(¬xi ∨ ¬xj) pairwiseO(N·K²)两组格互斥(¬ai ∨ ¬bi) for each value iO(K·pairs)2.3 直线组与相邻约束统一成“互斥对”蜂窝网格里最常见的约束是两类两条格子组成的“邻居对”互异以及若干格子组成的“直线组”互异。两者本质相同同一个值的两个出现不能同时落在已知格子里。所以共用一个函数void addDistinctPair(int a, int b) { for (int v 1; v K; v) { clauses.push_back({-litId(a, v), -litId(b, v)}); } } void addDistinctGroup(const vectorint group) { for (size_t i 0; i group.size(); i) for (size_t j i 1; j group.size(); j) addDistinctPair(group[i], group[j]); }a、b都是 cellId。每条子句的含义是格子 a 和格子 b 不能同时取同一个值 v。三条轴方向的直线组和相邻对在第 3 章生成后直接喂给这两个函数即可。到这里 CNF 主体已经齐了。接下来最容易被忽略的是初始提示数字——课程设计通常会给出若干已知数字这些不是复杂约束而是单文字子句clauses.push_back({litId(cell, given)});。如果愿意还可以顺手把该格其他候选值的否定也加进去让单位传播更快锁定解空间。3. 六边形网格的坐标与邻居表约束自动生成的正确性来源3.1 用轴向坐标生成 19 格蜂窝六边形网格最常见的存储方式是 axial coordinates两个坐标 q、r 就足够第三个 s 由s -q-r导出。六个邻居方向因此非常整齐不需要像二维数组那样做奇偶行判断。半径 R 的蜂窝形网格中心格到最远格距离 R总格数为3R(R1)1R2 时正好是 19 格也是这类题目最常用的规模。vectorCell cells; int R 2; for (int q -R; q R; q) { for (int r max(-R, -q - R); r min(R, -q R); r) { int s -q - r; cells.push_back({q, r, s, (int)cells.size()}); } } unordered_maplong long, int posToId; for (auto c : cells) posToId[(c.q 10) * 100 (c.r 10)] c.id;posToId的 key 是自己定义的编码q10和r10只是把负坐标平移成正数。实际项目里更干净的做法是mappairint,int,int或自定义哈希函数。生成后每个 Cell 的 id 从 0 排到 18后续所有约束生成都只谈 id不再碰坐标。3.2 按三个轴向提取直线组蜂窝数独里“行”的定义通常有三个方向q 相同、r 相同、s 相同。遍历这三个方向把相同坐标值的格子收集成一组长度大于等于 2 的组就是要施加“组内互异”的对象vectorvectorint groups; for (int axis 0; axis 3; axis) { for (int v -R; v R; v) { vectorint g; for (auto c : cells) { if ((axis 0 c.q v) || (axis 1 c.r v) || (axis 2 c.s v)) g.push_back(c.id); } if (g.size() 2) groups.push_back(g); } }这段代码会生成若干长度 5、4、3、2 的组。对每一组调addDistinctGroup就完成了“同一轴线数字不重复”的编码。调试时值得把组全部打印出来与手画网格对照R2 时每个轴方向取 -2..2 共 5 条线三条轴共 15 组组的总长度为 57。这个数字如果对不上说明坐标生成或分组条件有误不必等求解结果就能发现。3.3 邻居表生成边界处理决定正确性邻居关系用六个方向向量计算int dirQ[6] {1, -1, 0, 0, 1, -1}; int dirR[6] {0, 0, 1, -1, -1, 1}; for (auto c : cells) { for (int d 0; d 6; d) { int nq c.q dirQ[d]; int nr c.r dirR[d]; int ns -nq - nr; auto it posToId.find((nq 10) * 100 (nr 10)); if (it ! posToId.end()) { int nb it-second; if (c.id nb) // 每条边只加一次 addDistinctPair(c.id, nb); } } }c.id nb的判断避免同一条边被正反各加一次。越界坐标在posToId中查不到所以不会产生越界邻居。这个设计的正确性完全建立在 3.1 的坐标生成和 3.3 的 key 编码一致上只要 posToId 的 key 与查询 key 拼接方式不同邻居表会立刻缺边或多边。实验报告里我建议把邻居表打印出来人工抽查几格中心格应有 6 个邻居角格 3 个边格 4 个。这个检查能在五分钟内定位九成以上的建模错误。4. 只用C标准库实现DPLL单位传播、MRV决策与回溯4.1 子句存储与赋值栈不引入 MiniSat 等外部库时递归回溯实现一个可用 DPLL 完全可行。子句用vectorvectorint存每条子句内是正负文字。赋值状态开两个数组vectorint assign; // 0 未赋值1 为真-1 为假 vectorint reason; // 记录导致赋值的子句下标回溯时用于定位reason在简化版里不一定要实现完整回溯链但保留它对实验报告的“传播过程分析”很有用。回滚可以直接在递归进入前把整个 assign 拷贝一份子句数在千级别时拷贝开销可以忽略。这里的取舍是递归回溯追求正确性和可读性而不是和 MiniSat 比性能。4.2 单位传播传播函数反复扫描所有子句寻找“只有一个文字未满足”的子句把剩余文字赋成真。核心实现bool propagate(vectorvectorint clauses, vectorint assign) { bool changed; do { changed false; for (const auto cl : clauses) { int unassigned 0, lastPos -1; bool satisfied false; for (int lit : cl) { int val assign[abs(lit) - 1]; if (val 0) { unassigned; lastPos lit; } else if ((lit 0 val 1) || (lit 0 val -1)) { satisfied true; break; } } if (satisfied) continue; if (unassigned 0) return false; // 冲突 if (unassigned 1) { assign[abs(lastPos) - 1] (lastPos 0) ? 1 : -1; changed true; break; } } } while (changed); return true; }注意assign[abs(lit)-1]的下标litId 返回 1-based所以下标要减一。每次传播结束必须从头重新扫描因为新赋值可能立刻让前面的子句变成单位子句。这个 O(子句数 × 文字数) 的做法不求速度但正确性最直观。想提速可以课后实现 two-watched-literals那是另一个独立话题。4.3 决策策略MRV优先DPLL 的搜索效率主要由决策顺序决定。蜂窝数独规模小MRVMinimum Remaining Values足够每次传播后找一个候选值最少的格子再选它的候选值做分支。选格子的逻辑是int chooseVariable(const vectorint assign, int N) { int best -1, bestCnt K 1; for (int cid 0; cid N; cid) { int cnt 0; for (int v 1; v K; v) if (assign[varId(cid, v)] 0) cnt; if (cnt 1 cnt bestCnt) { bestCnt cnt; best cid; } } return best; }cnt 1是关键一个格子只剩一个可选值时传播阶段早就把它定下来了不需要决策。完整搜索函数是标准框架bool dpll(vectorvectorint clauses, vectorint assign) { if (!propagate(clauses, assign)) return false; int cid chooseVariable(assign, cells.size()); if (cid 0) return true; // 全部确定得到解 for (int v 1; v K; v) { if (assign[varId(cid, v)] ! 0) continue; auto saved assign; assign[varId(cid, v)] 1; if (dpll(clauses, assign)) return true; assign saved; } return false; }每次进入递归前拷贝 assign出了失败分支就整体恢复逻辑上最不容易出错。候选值的尝试顺序也可以做文章优先试“与已赋值邻居约束最多”的值能减少回溯次数。实验报告里可以对比有无该排序的运行数据作为性能分析素材。5. 跑通求解器输入格式、CNF检查与结果验证5.1 最小命令行实现主程序按三步走读题目、生成 CNF、求解并打印。题目文件用最朴素的文本格式每行格子编号 数字没有给出的格子不写3 5 7 2 18 9解出来之后把每个格子坐标和值一起打印方便对着图人工核对for (auto c : cells) { int val 0; for (int v 1; v K; v) if (assign[varId(c.id, v)] 1) { val v; break; } printf((%d,%d,%d) - %d\n, c.q, c.r, c.s, val); }搜索结束后 assign 数组保存的就是满足所有子句的赋值。这个输出格式还有个额外好处可以用脚本把 (q,r,s) 转成二维坐标画成六边形图直接在终端里看整个解的形状。5.2 用“不可满足”反向验证CNF正确性这是调试这类程序最有效的方法。选一个显然无解的题目比如相邻两个格子给了同一个数字或者同一个格子给了两个不同数字跑一遍。如果程序一秒内输出 UNSAT说明基本传播和回溯逻辑是通的如果输出了“解”CNF 生成或传播一定有 bug。再用空题目跑一遍应能快速输出一个平凡解。随后加 2~3 个提示数字做中等规模测试。一次典型运行结果如下测试输入变量数子句数决策次数结果无提示格171约 1700几十SAT相邻格给同值171约 17001UNSAT3 个提示格171约 1700几十SAT具体决策次数会随实现波动但数量级就是这样。实验报告放这种表比贴一大段运行日志更能说明程序行为。5.3 两类常见bug的排除方法第一类是坐标哈希冲突。(q10)*100 (r10)在半径小于 10 时没问题但如果把 R 改大或用更大网格这个编码可能碰撞。稳妥做法是用mappairint,int,int或者定义key(q,r)返回 pair。第二类是组约束漏了某个轴向的线。当某条线长度只有 1 时addDistinctGroup不会产生子句这是正确的但如果长度 2 的线因为 s 方向没收集到就会出现“同一行两个格子填相同数但求解器不报错”。遇到这种问题把 groups 全部打印出来逐条核对比看求解过程高效得多。6. 进阶把求解器改造成拓扑无关的通用工具6.1 用可配置的网格描述替代写死的19格目前的坐标生成把 R2 写死在代码里换一个网格形状就要改主流程。更通用的做法是把网络拓扑抽象成三张表格子列表、邻居对列表、分组列表。程序入口先读一个描述文件再交给同一套 CNF 生成和 DPLL// grid.txt N19 NEIGHBOR 0 1 NEIGHBOR 0 2 GROUP 0 1 2 3 4 GIVEN 3 5解析器读到NEIGHBOR就调用 addDistinctPair读到GROUP就调用 addDistinctGroup读GIVEN就生成单文字子句。这样同一个求解器可以解标准数独、锯齿数独、六角蜂窝。课程设计答辩时老师最常问“你这个能解别的题吗”这个改动只需要几十行代码就能从容回答。6.2 三个值得写进实验报告里的优化点第一把 pairwise 的“至多一”换成顺序编码sequential encoding每格子句数从 O(K²) 降到 O(K)第二用两个监视文字替代全表扫描传播让单位传播复杂度从 O(子句数) 降到 O(冲突传播链长度)第三冲突时做子句学习哪怕只记录最近一次冲突的原因子句UNSAT 判定也能快一个量级以上。前两个优化一晚上能写完第三个需要读几篇 SAT 入门材料但对 C 课程设计来说超出预期。注意别一上来就追求 MiniSat 级别性能19 格蜂窝数独的状态空间并没有大到必须靠工程优化才能收场。最终这套方案的价值在于把题目从“针对固定盘面的回溯程序”提升为“完整的 SAT 建模与求解演示”任意蜂窝形状、任意提示集合都能在笔记本上毫秒级判定可满足性。验证阶段把生成的 CNF 文件行数、各轴分组数、邻居边数打印出来和理论计算值对一遍比口头解释可靠得多。本文还有配套的精品资源点击获取