[C/C++] 舞蹈链算法实现：用C语言编写PostgreSQL UDF毫秒级求解数独

要实现在 PostgreSQL 中解决数独的“性能天花板”，必须编写一个基于C 语言的本地函数（UDF），并应用高德纳（Donald Knuth）提出的舞蹈链算法（Dancing Links, DLX）。

DLX 专门用于解决“精确覆盖问题”，它通过双向十字链表（四向指针）在搜索过程中动态地“切断”和“恢复”冲突节点，避免了数组拷贝或大规模内存移动，算法效率极高。

1. 编写 C 语言 UDF 代码 (sudoku_dlx.c)

这段代码实现了 DLX 的核心逻辑。它将 81 位数独转化为一个 729 x 324 的稀疏矩阵，并利用递归进行回溯搜索。

#include "postgres.h"
#include "fmgr.h"
#include "utils/builtins.h"
#include <stdbool.h>

PG_MODULE_MAGIC;

// 重命名为 DLXNode 避免与 PG 内核冲突
typedef struct DLXNode {
    int left, right, up, down;
    int col, row;
} DLXNode;

// 建议将大数组放在堆上或者作为静态变量，但在函数内初始化
static DLXNode dlx_matrix[100000];
static int col_count[325];
static int head, node_idx;
static int solution[81];

// 初始化 DLX 矩阵
static void init_dlx(int n) {
    for (int i = 0; i <= n; i++) {
        dlx_matrix[i].left = i - 1;
        dlx_matrix[i].right = i + 1;
        dlx_matrix[i].up = dlx_matrix[i].down = i;
        col_count[i] = 0;
    }
    dlx_matrix[0].left = n;
    dlx_matrix[n].right = 0;
    head = 0;
    node_idx = n + 1;
}

// 插入节点
static void add_node(int r, int c) {
    int first = -1;
    // 寻找当前行的第一个节点
    for (int i = node_idx - 1; i >= 325 && dlx_matrix[i].row == r; i--) {
        first = i;
    }
    dlx_matrix[node_idx].row = r;
    dlx_matrix[node_idx].col = c;
    dlx_matrix[node_idx].up = dlx_matrix[c].up;
    dlx_matrix[node_idx].down = c;
    dlx_matrix[dlx_matrix[c].up].down = node_idx;
    dlx_matrix[c].up = node_idx;

    if (first == -1) {
        dlx_matrix[node_idx].left = dlx_matrix[node_idx].right = node_idx;
    } else {
        dlx_matrix[node_idx].left = dlx_matrix[first].left;
        dlx_matrix[node_idx].right = first;
        dlx_matrix[dlx_matrix[first].left].right = node_idx;
        dlx_matrix[first].left = node_idx;
    }
    col_count[c]++;
    node_idx++;
}

static void cover(int c) {
    dlx_matrix[dlx_matrix[c].right].left = dlx_matrix[c].left;
    dlx_matrix[dlx_matrix[c].left].right = dlx_matrix[c].right;
    for (int i = dlx_matrix[c].down; i != c; i = dlx_matrix[i].down) {
        for (int j = dlx_matrix[i].right; j != i; j = dlx_matrix[j].right) {
            dlx_matrix[dlx_matrix[j].down].up = dlx_matrix[j].up;
            dlx_matrix[dlx_matrix[j].up].down = dlx_matrix[j].down;
            col_count[dlx_matrix[j].col]--;
        }
    }
}

static void uncover(int c) {
    for (int i = dlx_matrix[c].up; i != c; i = dlx_matrix[i].up) {
        for (int j = dlx_matrix[i].left; j != i; j = dlx_matrix[j].left) {
            col_count[dlx_matrix[j].col]++;
            dlx_matrix[dlx_matrix[j].down].up = j;
            dlx_matrix[dlx_matrix[j].up].down = j;
        }
    }
    dlx_matrix[dlx_matrix[c].right].left = c;
    dlx_matrix[dlx_matrix[c].left].right = c;
}

static bool search(int k) {
    if (dlx_matrix[head].right == head) return true;
    int c = dlx_matrix[head].right;
    for (int i = dlx_matrix[c].right; i != head; i = dlx_matrix[i].right) {
        if (col_count[i] < col_count[c]) c = i;
    }
    if (col_count[c] == 0) return false;
    cover(c);
    for (int i = dlx_matrix[c].down; i != c; i = dlx_matrix[i].down) {
        solution[k] = dlx_matrix[i].row;
        for (int j = dlx_matrix[i].right; j != i; j = dlx_matrix[j].right) cover(dlx_matrix[j].col);
        if (search(k + 1)) return true;
        for (int j = dlx_matrix[i].left; j != i; j = dlx_matrix[j].left) uncover(dlx_matrix[j].col);
    }
    uncover(c);
    return false;
}

PG_FUNCTION_INFO_V1(solve_sudoku_dlx);

Datum solve_sudoku_dlx(PG_FUNCTION_ARGS) {
    text *input_text = PG_GETARG_TEXT_PP(0);
    char *input = text_to_cstring(input_text);

    if (strlen(input) != 81) {
        ereport(ERROR, (errmsg("Input string must be exactly 81 characters")));
    }
    init_dlx(324);
    for (int r = 0; r < 9; r++) {
        for (int c = 0; c < 9; c++) {
            int num = (input[r * 9 + c] == '.') ? 0 : input[r * 9 + c] - '0';
            for (int v = 1; v <= 9; v++) {
                if (num == 0 || num == v) {
                    int row_id = (r * 81 + c * 9 + v);
                    add_node(row_id, r * 9 + c + 1);
                    add_node(row_id, 81 + r * 9 + v);
                    add_node(row_id, 162 + c * 9 + v);
                    add_node(row_id, 243 + ((r / 3) * 3 + c / 3) * 9 + v);
                }
            }
        }
    }
    if (search(0)) {
        char out[82];
        memset(out, 0, 82);
        for (int i = 0; i < 81; i++) {
            int r = (solution[i] - 1) / 81;
            int c = ((solution[i] - 1) % 81) / 9;
            int v = (solution[i] - 1) % 9 + 1;
            out[r * 9 + c] = v + '0';
        }
        PG_RETURN_TEXT_P(cstring_to_text(out));
    }

    PG_RETURN_NULL();
}

2. 编译与部署

要在 PostgreSQL 中使用此函数，可使用 Makefile：

MODULES = sudoku_dlx

PG_CONFIG = pg_config
PGXS := $(shell $(PG_CONFIG) --pgxs)
include $(PGXS)

确保 pg_config 命令在 PATH 环境变量中。

USE_PGXS=1 make
USE_PGXS=1 make install

以上代码已在 PolarDB for PostgreSQL 15 开源版本验证可用！

3. 创建 SQL 函数

在 SQL 终端中运行以下命令来定义函数，此操作在 数据库/中间件 的管理中很常见：

CREATE OR REPLACE FUNCTION solve_sudoku_dlx(text)
  RETURNS text
  AS 'sudoku_dlx', 'solve_sudoku_dlx'
  LANGUAGE C STRICT parallel safe;

4. 性能与效果对比

现在你可以直接在 SQL 中调用它了：

-- 测试“骨灰级”数独
SELECT solve_sudoku_dlx('8..........36......7..9.2...5...7.......457.....1...3...1....68..85...1..9....4..');

对比之前使用 PostgreSQL 内置的递归 CTE 解决同一个数独问题需要 700 多毫秒，现在使用 C语言 编写的 UDF 只需要约 1 毫秒。

postgres=# -- 测试“骨灰级”数独
SELECT solve_sudoku_dlx('8..........36......7..9.2...5...7.......457.....1...3...1....68..85...1..9....4..');
                               solve_sudoku_dlx
-------------------------------------------------------------------------------
 812753649943682175675491283154237896369845721287169534521974368438526917796318452
(1 row)

Time: 1.419 ms

使用正则表达式格式化输出更直观：

SELECT regexp_replace(
    board,
    '(.{9})',
    '\1' || chr(10),
    'g'
) from (values (solve_sudoku_dlx('8..........36......7..9.2...5...7.......457.....1...3...1....68..85...1..9....4..')))
 as t (board);

为什么这是最高效的？

1. 算法效率：DLX 算法的时间复杂度在数独问题上几乎是常数级的。即使是那些让递归 CTE 跑几百毫秒的题目，C 语言版的 DLX 通常在 1 毫秒（甚至更少） 内就能解完。
2. 内存开销：C 代码直接在栈和静态内存上操作，没有 PostgreSQL 递归查询时产生的庞大中间表（Temporary Tables）开销。
3. 无类型转换：位运算优化版 SQL 依然受限于 PG 的类型检查和函数调用开销，而 C 函数直接操作指针，消除了所有中间层。

性能进阶：SIMD优化探索

在 Apple M2 (ARM 架构) 上，针对舞蹈链算法（DLX）进行 SIMD (Single Instruction, Multiple Data) 优化是一个非常硬核的方向。

由于 DLX 算法的核心是基于稀疏矩阵的深度优先搜索（DFS），它涉及大量的链表指针跳转（Pointer Chasing）和递归回溯。这种非连续内存访问模式对 SIMD 并不友好。不过，通过改变思路，可以从以下几个维度压榨性能：

1. 内存布局优化

将 left, right, up, down 拆分为四个独立的结构化数组（SoA）。在频繁访问 up 和 down 的列覆盖操作中，可以提高 L1 Cache 的命中率。在初始化矩阵时，可以利用 NEON 指令集一次性填充多个索引数据。

2. 并行冲突检查

引入Bit-parallelism来加速。在进入 search 递归前，使用 NEON 指令并行预处理多个位置的候选数冲突检查。

3. 利用硬件加速

在寻找 col_count 最小的列时（启发式搜索），利用 AArch64 的 CNT 指令或 __builtin_popcount 硬件加速统计操作。

4. 消除递归开销

将递归改写为循环，并使用一个紧凑的固定数组手动管理栈，减少函数调用和分支预测压力。

5. 代码示例（NEON预处理）

可以在函数开头加入 NEON 指令进行快速初始化：

#include <arm_neon.h>

// 示例：使用 NEON 指令并行清零列计数器
void fast_zero_columns(int *counts) {
    uint32x4_t zero = vdupq_n_u32(0);
    for (int i = 0; i < 324; i += 4) {
        vst1q_u32((uint32_t *)&counts[i], zero);
    }
}

6. 系统级优化

• 使用 posix_memalign 确保内存分配与 64 字节缓存行对齐。
• 手动展开 DLX 矩阵中每行固定 4 个节点的循环。

结论与适用场景

在 M2 上，原版 C 实现的 DLX 解决“骨灰级”数独通常耗时在 30-100 微秒 左右。

• SIMD 收益：可能将耗时进一步压低 20%-30%，但会显著降低代码可读性。
• 更优策略：对于每秒处理数万个不同数独题目的场景，使用多核并行（如OpenMP）或多连接并发的收益远大于复杂的 SIMD 优化。

DBA 提示

• 安全性：在生产环境加载 C 语言动态库需要超级用户权限，因为它直接运行在数据库进程的内存空间中。请确保代码经过严格测试，以防段错误（Segment Fault）导致数据库服务中断。
• 适用场景：这种方式最适合需要处理海量数独数据（如每天几百万道题）的高性能计算场景。