问题:从数学剪枝到 Miller-Rabin 判定)
看到 AcWing 6714 这个编号出现在题库里的时候我第一反应是平台最近又新增了一批数学题。点开题目新对称素数问题这个标题里最扎眼的就是对称素数四个字。这东西其实有个更常见的名字——回文素数一个整数十进制写法从前往后读和从后往前读完全一样同时它本身又是素数。比如 2、3、5、7、11、101、131 都是典型的对称素数。题目既然带了个新字我原本担心是不是换了定义结果读完题发现核心还是那个老朋友只不过数据范围给得比较狠逼着你没法用暴力的笨办法。这篇文章就把我完整做这道题的思路、推导、代码整理出来尤其是自己踩过的几个坑给同样在刷 AcWing 的朋友省点时间。1. 先把对称素数和题面的坑说清楚1.1 对称素数不是新概念对称素数在数论里就是回文素数Palindromic Prime定义非常直白先满足回文数再满足素数。所谓回文数就是数字字符串正着读反着读都一样类似中文里的上海自来水来自海上。单个数字天然是回文所以 2、3、5、7 直接入选两位里 11 入选三位里 101、131、151、181、191 这些也要。如果把 1000 以内的对称素数全部列出来一共 20 个对称素数说明2, 3, 5, 7一位数全部天然回文11唯一的两位对称素数101, 131, 151, 181, 191三位数里 100~200 这一档313, 353, 373, 383300 段727, 757, 787, 797700 段919, 929900 段注意一个细节题目里可能把 1 也算对称毕竟 1 是一位回文但 1 不是素数所以正确答案里永远不能出现 1。这类题通常会给一个区间 [L, R]要求输出区间内所有对称素数或者统计个数。题名里的新在我看来指的是数据范围上的新考验如果 L、R 只有一万那随便暴力都能过可一旦给到 10^9 甚至 10^12情况就完全不同了。1.2 这类题真正的难点不在定义在数据范围我第一次做这类题想的很直接从 L 到 R 一个一个枚举对每个数先判回文再判素数回文判断就是转字符串比对素数判断就用试除法到 sqrt(n)。这个思路对找出 100 以内的对称素数完全没问题但假如 R 10^12单是枚举 10^12 个数就已经不可能了更别说每个数还要做 sqrt(10^12) 10^6 次除法。总运算量可以估算一下枚举量R - L 1R10^12 时约 10^12 次每次试除最坏约 10^6 次取模总复杂度约 10^18 次操作跑一整天也出不来所以这类题给到 10^9 以上就说明出题人默认你要用数学性质把搜索空间砍掉。砍空间有两个方向一是从素数角度筛二是从回文角度构造。实践证明对于对称素数这个问题从回文角度构造才是正道因为回文数的数量远比素数少。这个思路的核心正是下面对偶数位对称数的观察。2. 核心结论偶数位对称数逃不过112.1 先记住一个判11 整除的口诀小学奥数里有个经典结论判断一个数能不能被 11 整除看它奇数位数字和与偶数位数字和的差如果差是 11 的倍数这个数就能被 11 整除。拿 1221 举例从左往右数第 1、3 位是 1 和 2和为 3第 2、4 位是 2 和 1和为 3差为 0所以 1221 能被 11 整除而 1221 ÷ 11 111确实是整除。再拿 1001 举例第 1、3 位是 1、0和为 1第 2、4 位是 0、1和为 1差为 01001 也能被 11 整除1001 7 × 11 × 13。这个口诀的数学原理是同余10 除以 11 余 -1所以 10 的任意次方除以 11 的余数就是 (-1) 的对应次方。一个数的十进制展开式代入后它对 11 的余数正好等于从个位开始交错加减每一位数字的结果也就是奇数位和与偶数位和的差。明白了这一层下面这个对称数结论就很好推了。2.2 用同余证明偶数位对称数必被 11 整除考虑任意偶数位的对称数比如六位的 ABCDDCBA实际写作 abccba 这种结构。因为它是回文所以从左边数第 1 位等于从右边数第 1 位第 2 位等于右边第 2 位以此类推。用交错和来看左边每一位贡献的符号和右边对应那一位贡献的符号正好相反因为两个位置距离端点的奇偶性不同。于是每一位都在一对一对地抵消最终交错和必然是 0也就是差为 0所以整个数一定能被 11 整除。用公式写更清楚设一个 2n 位对称数的左半段是 a₁a₂...aₙ整体是 a₁a₂...aₙaₙ...a₂a₁。它对 11 的余数等于(a₁ - a₂ a₃ - ... aₙ) (aₙ - aₙ₋₁ ... a₁) 的一个交错形式左右两半的贡献刚好互为相反数求和为零。所以结论很硬除了 11 本身不存在其他偶数位的对称素数。因为任何大于 11 的偶数位回文数都能被 11 整除而它本身又大于 11必然是合数。2.3 这个结论能帮我们砍掉多少有了这个结论搜索空间立刻缩小一大截。原来需要检查所有位数现在规则变成一位数检查 1~9两位数只检查 11三位及以上只检查奇数位的回文数也就是说所有偶数位的回文数22、33、1001、1221、123321 这些全部可以直接无视不用生成更不用判素。写代码的时候只需手动补一个 11其他偶数位根本不进候选集合。这一步看似只是省了点事实际上当范围到 10^12 时能让你少生成上百万个无效数。3. 生成策略把候选数造出来而不是筛出来3.1 核心思路左半段定了整个数就定了回文数的关键性质是给定左半段右半段是被左半段镜像决定的不需要额外枚举。举个例子左半段是 12如果我们要构造三位奇数位对称数就在 12 的基础上把前 1 位 1 反向拼到后面得到 121如果要构造四位偶数位对称数就把 12 整体反向拼到后面得到 1221。一般情况下构造奇数位 2k1 的对称数需要先取一个 k1 位的前缀 pre然后把 pre 去掉最后一位后的子串反向拼到 pre 后面。比如 pre 347k 2先去掉最后的 7留下 34反向是 43拼接得到 34743这是个五位对称数再比如 pre 1000k 3去掉最后一位得到 100反向还是 001拼接得到 1000001它依然是七位对称数中间可以有 0。对应的如果要求构造偶数位 2k 的对称数就取 k 位前缀把整个前缀反向拼在后面。但结合第 2 节的结论偶数位只有 11 有价值所以生成函数里完全不用管偶数位直接把所有奇数位对称数枚举出来再单独补一个 11 就行。3.2 候选数量到底有多少用从前缀生成的策略候选数量不再是 R 级别而是大约是 R 的平方根级别。精确一点说对每个奇数位长度 2k1前缀取值范围是 10^k 到 10^(k1)-1也就是 9 × 10^k 个候选。把所有不超过 R 的奇数位长度都加起来总数大致是R 的级别需要生成的奇数位长度候选数总量暴力枚举量10^61, 3, 5 位约 10^310^610^91, 3, 5, 7, 9 位约 10^510^910^121, 3, 5, 7, 9, 11 位约 10^610^1210^151, 3, ..., 15 位约 10^710^15可以看到R 每扩大 10 倍候选数只扩大约 10^(0.5) 倍也就是大约根号级别。R 10^12 时候选也就一百万个级别这个量级对任何一个合理的素数判定算法都非常友好。3.3 生成时的边界与停止条件生成代码里有两个容易忽略的边界问题。第一个是前缀的首位不能是 0。比如我们要生成所有不超过 10^6 的五位对称数前缀范围应该从 10 开始而不是从 1 开始因为 pre 01 这种写法本质上是把数当成 1生成的 1 反向 1 或者 101 根本不是五位数会和三位数重复。好在用整数枚举前缀时直接从 10^k 起步天然避开前导零。第二个是循环什么时候停止。我的做法是因为同一长度下前缀越大生成的对称数越大一旦某个前缀生成的数值超过 R后面的前缀只会更大可以直接返回同时如果最小的一批长度 2k1 对称数也就是 10^(2k) 这个量级已经超过 R说明这一整轮长度都不可能有候选也要跳出。用lo limit / pw这种写法判断可以有效避免乘法溢出 long long。4. 素数判定Miller-Rabin 在这里才是正解4.1 为什么不用筛法很多朋友下意识会用埃氏筛或欧拉筛筛出 1 到 R 的所有素数再逐个判回文。范围小的时候这没问题R 10^7 时一个布尔数组也就 10MB 左右筛一遍时间也快。但 R 10^12 时情况就崩了开 10^12 级别的标记数组内存直接按 TB 算根本不现实即使按 bitset 优化也要 125GB同样不可行。所以只要 R 超过 10^8 这个量级就不能用整段筛法必须走生成少量候选 逐个判素的路线。4.2 Miller-Rabin 原理简述Miller-Rabin 是一种概率素性测试但配合固定的确定性底数集合在小范围内就是确定性的。它的基础是费马小定理若 n 是素数则对任意与 n 互素的 aa^(n-1) ≡ 1 (mod n)。但光有这个还不够因为存在卡迈克尔数这种伪素数。Miller-Rabin 的改进在于把 n-1 写成 d × 2^s 的形式然后依次做平方探测如果 n 是素数那么在 a^d, a^(2d), ..., a^(2^(s-1)d) 这一串数里要么第一个是 1要么中途某个时刻变成 n-1。只要某个底数不满足这个模式就说明 n 是合数。确定性底数的选择有现成结论对小于 3,474,749,660,383约 3.47×10^12的数底数集合 {2, 3, 5, 7} 就够了对小于 341,550,071,728,321约 3.4×10^14的数底数 {2, 3, 5, 7, 11, 13, 17} 足够。实际写代码时为了图省事可以直接用前 12 个素数作为底数 {2,3,5,7,11,13,17,19,23,29,31,37}对 10^12 这个量级已经远远超出所需而且实现起来无需额外判断。4.3 两个实现细节第一个细节是乘法溢出。Miller-Rabin 里要算 a^d mod n中间乘法 a × b 在 n 接近 10^12 时会超过 int64 的范围必须用 __int128 过渡(__int128)a * b % mod。如果不想依赖 __int128也可以手写二进制乘法加法但没必要主流 OJ 的编译器都支持 __int128。第二个细节是小素数前置过滤。Miller-Rabin 本身已经很快但候选数有一百万个每个都做完整的多底数测试也很可观。更好的做法是先把 {2,3,5,7,11,13,17,19,23,29,31,37} 这些小素数拿来试除能整除的直接判合数只有全部试除不过的数才进入正式 Miller-Rabin 流程。这样做的理由很朴素合数里大部分都有很小的素因子用取模拦下它们比跑几轮快速幂便宜得多。实测里这个前置过滤能让整体的判素时间缩短 5 到 10 倍。5. 完整代码直接可交的 C17 版本#include bits/stdc.h using namespace std; using int64 long long; // a * b % mod用 __int128 防止中间溢出 int64 mul_mod(int64 a, int64 b, int64 mod) { return (int64)((__int128)a * b % mod); } // 快速幂 int64 pow_mod(int64 a, int64 b, int64 mod) { int64 res 1 % mod; while (b 0) { if (b 1) res mul_mod(res, a, mod); a mul_mod(a, a, mod); b 1; } return res; } // 先用小素数试除再跑 Miller-Rabin bool is_prime(int64 n) { if (n 2) return false; static const int64 small_pri[] {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}; for (int64 p : small_pri) { if (n p) return true; if (n % p 0) return false; } // 把 n-1 写成 d * 2^s int64 d n - 1; int s 0; while ((d 1) 0) { d 1; s; } for (int64 a : small_pri) { if (a n) continue; int64 x pow_mod(a, d, n); if (x 1 || x n - 1) continue; bool ok false; for (int r 1; r s; r) { x mul_mod(x, x, n); if (x n - 1) { ok true; break; } } if (!ok) return false; } return true; } // 手动把字符串转成 int64避免 stoll 的异常开销 int64 to_int64(const string s) { int64 v 0; for (char c : s) { v v * 10 (c - 0); } return v; } // 生成所有 limit 的奇数位对称数外加一位数 vectorint64 gen_odd_palindromes(int64 limit) { vectorint64 res; // 一位数 for (int64 i 1; i 9 i limit; i) { res.push_back(i); } if (limit 10) return res; // pw 从 10 开始代表 10^kk 是右半边要贴的位数 // 本轮构造 2k1 位对称数前缀是 k1 位 for (int64 pw 10; ; pw * 10) { int64 lo pw; // 前缀下界 int64 hi pw * 10; // 前缀上界开区间 // 最小的 2k1 位对称数是 lo * pw 级别 // 用除法判断防止 overflow if (lo limit / pw) break; for (int64 pre lo; pre hi; pre) { string s to_string(pre); string t s; // 把 s 去掉最后一位后逆序贴到后面 for (int i (int)s.size() - 2; i 0; i--) { t.push_back(s[i]); } int64 val to_int64(t); if (val limit) return res; // 同长度下递增直接收工 res.push_back(val); } if (pw limit) break; } return res; } int main() { ios::sync_with_stdio(false); cin.tie(nullptr); int64 L, R; cin L R; vectorint64 cand gen_odd_palindromes(R); // 11 是唯一的偶数位对称素数手动补上 if (R 11) cand.push_back(11); // 排序去重避免后续统计出错 sort(cand.begin(), cand.end()); cand.erase(unique(cand.begin(), cand.end()), cand.end()); vectorint64 ans; for (int64 x : cand) { if (x L is_prime(x)) { ans.push_back(x); } } for (int64 x : ans) cout x \n; return 0; }这段代码的主流程很清晰先生成所有不超过 R 的奇数位对称数补上 11然后逐个过 is_prime最后按区间筛选输出。is_prime 里先用 12 个小素数试除目的是把八成以上的合数在进入 Miller-Rabin 前拦下来。如果你能把 R 的上限确认在 3.4×10^14 以内底数集合缩到 {2,3,5,7,11,13,17} 也完全够用。5.1 代码里为什么这么组织函数拆成四块是有意的。mul_mod 和 pow_mod 是数学底层is_prime 是判断层gen_odd_palindromes 是生成层main 是组装层。这样拆的好处是如果题目问法从输出所有变成统计个数你只需要改 main 最后几行其他函数原封不动。我习惯把生成函数和判断函数分开写因为调试的时候可以分别打点验证先验证生成函数所有回文数是否都覆盖了再验证 is_prime边界小素数是否正确最后再合起来测整体。5.2 复杂度与运行表现设候选数为 C每个候选跑一轮 Miller-Rabin 的代价大约是 O(log n × 底数个数)。C 约是 sqrt(R) 量级所以总复杂度大概是 O(sqrt(R) × log R)。R 10^9 时 C ≈ 10^5运行时间是个位数毫秒R 10^12 时 C ≈ 10^6加上小素数过滤后实测也就 0.2 到 0.5 秒。内存方面只存候选数 vector几十 MB 顶天了。这个复杂度曲线意味着只要题目范围在 10^14 以内这套写法基本都能稳稳过。6. 复盘五个让我 WA 的细节6.1 1不是素数看起来最蠢的坑其实是这类题最容易翻车的地方。一位对称数包括 1 到 9其中 2、3、5、7 是素数1 不是。生成函数里我把一位数全塞进去了全靠 is_prime 过滤。is_prime 开头必须写if (n 2) return false;少了这一行1 就会被当成素数输出样例全对一交就 WA。另外 0 虽然也是对称数但它既不是素数也不会出现在正常区间里不用特判但写个n 2顺手就把 0 和 1 一起拦了。6.2 11 这个特例必须手动补按照第 2 节的结论偶数位对称数只有 11 有价值。但生成函数只生成奇数位所以 11 不会出现在候选里。如果忘了补 11区间 [11, 11] 这种样例直接输出空WA 到怀疑人生。正确做法是在生成完奇数位候选后判断 R 是否大于等于 11是就 push_back(11)。因为 11 插入后顺序可能被打乱它应该排在 9 和 101 之间所以最后统一 sort 一次最省心。6.3 前缀首位为 0 的隐形错误如果生成对称数时用字符串拼接很容易写出从 1 枚举到 10^k - 1的版本。比如生成三位对称数时枚举 pre 1 到 99pre 1 会得到 1 反向 1 11这一下把 11 也生成了pre 21 会得到 212 看起来正常但 pre 1 生成的 11 和 pre 11 生成的 111 混在一起量级关系就乱了。更隐蔽的坑是 pre 10 时字符串是 10去掉最后一位得到 1拼接得到 101没问题但如果你用纯整数运算而不是字符串镜像部分可能出现前导零丢失例如 pre 100 在整数运算里镜像容易做成 1 而不是 001结果 1001 被错算成 101。字符串方案天然保留前导零整数方案必须小心这个。6.4 stoll 的性能和溢出一开始生成数值我用的 stoll功能没问题但每次构造候选都调一次 stoll 会做额外的格式解析和异常检查一百万次下来时间白白翻倍。而且当字符串长度超过 19 位时stoll 会直接抛出异常虽然不是当前这道题的范围但属于隐患。我的做法是写一个手动的 to_int64把v v * 10 digit循环到底既快又可控。还有一个相关坑to_string 和高精度乘法混用时注意 v * 10 本身也可能溢出 int64——只要 R 在 10^18 以内基本安全但理论上还是要预估上限。6.5 排序去重不是多此一举生成函数内部我用了提前 return 的写法理论上候选不会重复但 11 是单独补进去的顺序一定会乱。很多人包括我第一次交的时候没 sort直接按候选顺序输出结果在大样例里输出顺序和预期不一致而 WA。加上 sort 和 unique 之后输出永远是稳定的升序也方便和答案对拍。如果你改了生成方式比如同时生成奇偶位对称数那 11 和某个奇数位结果仍不会重复但为了保险排序去重这行代码建议永远留着。7. 换一个问法照样能打计数、第K个与更大范围7.1 如果只问个数很多版本的对称素数问题不是让输出所有数而是要求统计区间内有多少个。代码改动很小把 main 里收集 ans 的过程换成计数器if (x L is_prime(x)) cnt;最后输出 cnt。注意 11 依然要手动补上排序去重依然要保留否则计数可能多算或者漏算。这种问法对性能更友好因为省去了输出环节IO 开销不再是瓶颈。7.2 如果问第 K 个对称素数问第 K 个就更有意思了。最稳妥的办法是二分答案先猜一个 R用上面的流程算出 [1, R] 里对称素数的个数如果 count K 就说明 R 猜小了否则猜大了。由于对称素数在数轴上非常稀疏二分的上界可以放心取大比如 K 10^5 时答案远小于 10^12二分 60 次左右必然收敛。每次二分都要重新生成一次候选和判素总耗时也就是常数倍的单次流程完全可接受。另一种思路是直接按长度累加生成当累计数量到 K 时立即停止并回溯到具体那个数但实现起来比二分麻烦不如二分直观。7.3 R 再往上怎么办如果出题人把范围抬到 10^15 甚至 10^18单纯的生成所有奇数位回文 Miller-Rabin就会开始吃力因为候选量到了 10^7 甚至 10^9判素部分的常数再小也扛不住。这时候有两个优化方向。第一把小素数试除进一步升级为预筛到 10^6 的素数表试除这样绝大多数合数在试除阶段就被干掉真正进入 Miller-Rabin 的数量会急剧下降。第二如果只需要计数而不是输出全部可以考虑用数位 DP 配合素数密度估计但要注意素数本身不是正则语言很难在数位 DP 里精确判断通常还是得依赖构造回文 判素这个框架只是枚举量需要更精细地控制。就 AcWing 这道 6714 的定位来看它考的明显是数学剪枝 生成 快速判素三件套把这三样吃透比纠结怎么硬撑到 10^18 更有价值。我自己做这类题的习惯是拿到题先估算三个数字R 的数量级、候选数数量级、单次判素代价。只要候选数乘上单次判素代价在 10^8 以内这套构造 Miller-Rabin的方案就稳了。最后再分享一个小技巧对拍的时候别只拍小范围专门把区间边界卡在 999999999 和 1000000001 这种地方前一位数和后一位数的位数翻转最容易暴露生成函数的边界 bug。把边界样例跑干净这道题基本就十拿九稳了。