把 USTC 计算数论这门课的笔记、实验代码和错题本重新翻出来整理了一遍前后断断续续花了两周时间。越整理越觉得这门课值得单独写一篇记录因为它跟本科阶段那种证明存在性的数论完全不是一条路子。初等数论考的是你能不能用反证法、能不能构造出一个优雅的同余式计算数论考的是你能不能在可接受的时间、可接受的内存里把一个具体的、动辄几百位的大整数问题真正算出来。素性检测、整数分解、离散对数、格约化、椭圆曲线点计数这些词听起来都很密码学但落到纸面上全是可以写成代码、可以跑出具体数字、可以用复杂度公式标出代价的算法工程。这篇记录适合三类人正在上这门课、被实验报告折磨的同学想自学 USTC 计算数论但不知道从哪个点切入的人以及数论早就忘光、只想看看这些算法到底怎么跑起来的工程师。我会把课程的整体框架、每个算法我自己理解它的方式、能直接抄走的代码、以及踩过的坑全部写出来。代码都在 Python 配 gmpy2/SageMath 的环境下验证过参数选择也都给了理由不是随手抄的教科书结论。1. 课程整体设计与我的学习路线拆解1.1 计算数论到底在算什么可以先给这门课画一个很朴素的边界它研究的是数论对象的可计算性、算法与复杂度。同样是问这个数是不是素数初等数论会用欧拉定理去证伪计算数论会直接估算给定 2^1024 量级的输入你要付出多少位运算。所以整门课的视角其实是计算机的视角输入的规模不是数值本身而是它的二进制长度 log₂n一切复杂度都按位运算来算。这个转换我刚上第一周的时候特别不适应总觉得分解 80 位的数能有多难直到真的动手去跑才发现指数级和亚指数级之间的差距有多残酷。课程内容大致压在三个核心问题上判定素性、整数分解、离散对数。这三个问题构成了整个计算数论的骨架剩下的有限域算术、椭圆曲线、格基约化本质上都是围绕它们延伸出来的工具与变体。有意思的是这三个问题的难度并不对称素性检测已经是多项式时间的确定性算法AKS整数分解至今只有亚指数时间的通用算法离散对数在有限域上的难度和分解问题有很深的关联但椭圆曲线上的离散对数目前没有已知的亚指数算法。这种难度不对称本身就是课程里反复强调的东西。我个人的感受是这门课真正的门槛不在数学而在工程感。你要能判断一个算法在什么规模下会崩要能根据输入的样子选算法——比如待分解的数如果 p-1 光滑那先上 Pollard p-1 比直接上 rho 快得多。这种看菜下饭的判断力是课本给不了的只能靠实验一个个跑出来。1.2 课程内容地图五个模块的递进关系我把整门课的内容重新梳理成五个模块这个划分不是课表上的是我自己复习时整理出来的但非常好用第一个模块是基础算术层。涵盖大整数表示、加减乘除、模运算、模幂、扩展欧几里得、中国剩余定理、连分数与有理逼近。听起来最无聊实际上是后面所有算法的性能地基。模幂写错了后面全部白搭。第二个模块是素性与分解。从试除法、费马小定理、Carmichael 数讲起然后是 Miller-Rabin、Solovay-Strassen、Baillie-PSW再往下是 Pollard rho、Pollard p-1、ECM椭圆曲线分解法、二次筛法最后简单介绍数域筛法的思想。第三个模块是离散对数。Baby-step giant-step、Pohlig-Hellman 归约、Index Calculus以及椭圆曲线上的 BSGS 变体。第四个模块是代数结构上的计算。有限域的构造与运算、Legendre/Jacobi 符号、二次剩余与 Tonelli-Shanks 求模平方根、原根与阶的计算。这部分是工具库很多题目都是靠它们化简的。第五个模块是格与二次型。格基约化LLL、短向量问题、Coppersmith 方法的思想。这部分是课程后半段的分水岭也是很多人掉队的地方因为思路跟前四块完全不同——前四块是数的运算第五块是向量空间的运算。这五个模块的递进是很清楚的一提供工具二是核心问题三是对偶问题四是把问题搬到更一般的代数结构上五是把问题搬到几何/线性代数的语言里。复习的时候如果按这个顺序串比按章节顺序啃要顺得多。1.3 数学前置门槛与我的补课顺序课程对前置知识的要求是名义上不高实际上不低。名义上只需要初等数论加一门编程语言但实际上手写代码时会发现抽象代数的缺口会立刻暴露出来有限域 GF(p^k) 的构造需要你理解多项式环和理想椭圆曲线上的点加法需要你理解群、域扩张、射影坐标LLL 的正交化用的是 Gram-Schmidt。这些东西如果不熟代码写出来是对的但你不知道它为什么对一旦出问题就完全没法调试。我的补课顺序是这样的供参考先把初等数论的四个定理吃透费马小定理、欧拉定理、威尔逊定理、中国剩余定理这几条是所有算法的公共前提然后花三天补抽象代数里跟本课程相关的部分重点是群的定义与阶、有限域的结构、多项式环最后再补一点概率论因为 Miller-Rabin 和 Pollard rho 都是概率算法你需要理解误判概率是什么意思以及为什么重复 k 次能把误判概率压到 4^(-k)。编程这块我的建议是直接用 Python 起步不要一上来就写 C 或者 Rust。原因很实际这门课的实验重点是算法逻辑和复杂度验证不是内存管理和位操作。Python 的 int 天然就是任意精度pow(a, b, m)自带快速模幂你把精力全花在算法本身就好。等到你的感觉对了、要做性能对比了再换 gmpy2 或者直接上 C那才是有意义的优化。2. 核心算法细节解析与实操要点2.1 大整数运算所有性能问题的地基很多人低估了底层算术的重要性觉得乘法就是乘法。但在计算数论里乘法的复杂度直接决定了上层算法的可行规模。朴素的小学乘法是 O(n²)Karatsuba 是 O(n^1.585)Toom-Cook 三路大概是 O(n^1.465)再往上用 FFT 可以做 O(n log n log log n)。为什么这很重要因为 Pollard rho 里每一次迭代都要做一次模乘二次筛法里要做几十万次模乘底层的常数差异会被上层放大成能跑出来和跑不出来的区别。Python 内置 int 在几百位这个量级已经够用底层对中等规模的小整数其实用的是 Karatsuba 一类的算法大整数才会走更高级的乘法。gmpy2 的优势不只在于快更在于它提供了很多数论原语gcd、invert、powmod、is_prime、next_prime。做实验的时候如果你用纯 Python 手写这些代码量会大很多而且容易在边界情况上出错。这里有个我踩过的坑必须提醒别用浮点数判断大整数。math.isqrt是安全的但如果代码里出现int(n ** 0.5)当 n 超过 2^53 之后精度就不可靠了会造成试除范围不够、判平方失败这类幽灵 bug。凡是需要开方的地方一律用gmpy2.isqrt或math.isqrt。2.2 素性检测从试除到 Miller-Rabin素性检测的演进路线特别能体现计算数论的思维方式。试除法是 O(√n)对 10^20 量级完全不可行费马小定理说如果 n 是素数那么对所有与 n 互素的 a 都有 a^(n-1) ≡ 1 (mod n)问题是它的逆命题不成立——Carmichael 数比如 561 3×11×17对所有与它互素的底 a 都能通过费马测试这类数虽然稀疏但无穷多所以费马测试不能当素性判据用。Miller-Rabin 的关键改进是把 n-1 写成 2^s·dd 是奇数然后考察 a^d, a^(2d), a^(4d), ..., a^(2^(s-1)·d) 这一串值。如果 n 是素数这一串里要么第一项就是 1要么在某个位置出现 -1即 n-1。如果都不满足n 一定是合数如果满足n 只是可能是素数。单次测试把合数误判为素数的概率不超过 1/4所以重复 k 次、每次独立随机取底误判概率就压到 4^(-k)。取 k20误判概率大约是 10^(-12)工程上完全够用。但更漂亮的结论是确定性版本如果 n 小于 3,317,044,064,679,887,385,961,981约 3.3×10^24那么用前 13 个素数 {2,3,5,7,11,13,17,19,23,29,31,37,41} 作为底Miller-Rabin 的判定是确定性的没有任何误判。这个界在工程里极其实用——绝大多数场景下的整数都小于它你可以直接用固定基既快又绝对正确不需要随机数也不需要概率论证。我在实验报告里专门验证过这个界附近的几个数确实是严丝合缝的。还有几个值得一提的替代方案Solovay-Strassen 用 Jacobi 符号判定理论上更早但实际比 Miller-Rabin 慢课程里主要作为历史线索讲AKS 是第一个多项式时间的确定性素性检测算法理论意义巨大但常数大得离谱实际完全不能用Baillie-PSWMiller-Rabin 底 2 强 Lucas 测试在实践中从没被发现有反例工程上常被当作事实上的确定性算法。注意做实验的时候不要用random.random()取底那样拿到的是浮点数。必须用random.randrange(2, n-1)而且要确认 n-1 这个上界本身没问题。2.3 整数分解Pollard rho 与 p-1 的选型逻辑分解是整门课最耗时间的部分也是最讲究策略的部分。先建立一个基本认知没有任何已知的多项式时间通用分解算法。目前最好的通用算法是数域筛法复杂度是亚指数级的 L_n[1/3, 1.923...]这意味着你可以分解一两百位的数但再往上就要付出巨大的算力代价。实际做题时用的都是针对特定结构的算法核心思路是根据 n 的因子可能长什么样来选刀Pollard p-1 针对的是n 有一个素因子 p且 p-1 是 B-光滑的。它的原理是这样的如果 p-1 的所有素因子都不超过 B那么 p-1 整除 B!或者整除 lcm(1..B)于是对任意与 p 互素的 a有 a^(B!) ≡ 1 (mod p)也就是说 p 整除 gcd(a^(B!) - 1, n)。只要这个 gcd 不等于 n 本身就成功拿到一个非平凡因子。这个算法最大的优点是快只要 p-1 光滑几秒钟就能出结果缺点是 p-1 只要有一个大素因子就彻底失效。Pollard rho 是无结构依赖的通用启发式算法。它的思路非常巧妙用 f(x) x² c mod n 构造一个伪随机序列因为模 n 的序列会碰撞但真正有用的是它在模 p 意义下的碰撞p 是 n 的未知因子。由生日悖论模 p 的序列大概在 O(√p) 步之后就会碰撞。用 Floyd 判圈法同时跑 x 和 f(f(x))每次算 gcd(|x - y|, n)一旦这个 gcd 大于 1 就抓到了一个因子。期望复杂度是 O(n^(1/4))也就是对一个 2^100 的数大概 2^25 次迭代这在现代机器上也就几秒。选型上的经验是先跑 p-1成本极低试一下不亏再跑 rhorho 失败就换 c 重跑实在不行再考虑二次筛法。还有一个很实用的加速技巧是把若干次 gcd 合并起来算——连续迭代 100 次才做一次 gcd把这 100 个差值的乘积一起拿去做 gcd这样能省掉大量昂贵的 gcd 调用。这个技巧在 Brent 的改进版本里被发挥到了极致。2.4 离散对数BSGS 与 Pohlig-Hellman 的配合离散对数问题问的是给定群 G、生成元 g 和元素 h求 x 使得 g^x h。在有限域乘法群和椭圆曲线群上这都是困难问题也是很多密码方案的安全基础。课程里讲的方法按群的大小分两档。Baby-step giant-step 是通用的 O(√n) 时空算法思路是中间相遇。令 m ⌈√n⌉把 x 写成 x im j0 ≤ j m那么 g^(imj) h 等价于 g^j h·(g^(-m))^i。于是先把所有 g^j小步存进哈希表再枚举 i 去查表大步。一旦匹配上就得到 x。它的问题是要存 √n 个元素n 大到一定程度内存就爆了。Pohlig-Hellman 解决的是群阶是光滑数的情况。如果群的阶 N 分解成 ∏ p_i^(e_i)那么可以先把离散对数问题分别归约到每个素数幂阶的子群上求解再用中国剩余定理拼回去。这样一来总的复杂度从 O(√N) 降到 O(∑ e_i(log N √p_i))。在阶光滑的时候这个提升是数量级的。所以标准的做法是先用 Pohlig-Hellman 把阶分解掉剩下的大素数阶部分才用 BSGS。在有限域上还有 Index Calculus 这类亚指数算法思路类比分解里的筛法先选一个因子基收集关系式最后解线性方程组。椭圆曲线上目前没有对应的亚指数算法这也是椭圆曲线在同等安全强度下密钥更短的原因。2.5 二次筛法与格约化课程的分水岭二次筛法是第一个真正能分解大数的算法思想是寻找 x² ≡ y² (mod n) 且 x 不等于 ±y (mod n)此时 gcd(x - y, n) 就会给出非平凡因子。具体做法是选一个因子基一堆小素数把 ⌈√n⌉ i 的平方对 n 取模看结果能不能在因子基上完全分解光滑数能分解的就收集起来每一行记录一个 0/1 向量表示各素因子的奇偶性最后在 GF(2) 上找这些向量的线性相关组合。一旦找到一组相关组合就得到了 x² ≡ y²。这里最难的其实不是筛法本身而是线性代数。收集到几万到几十万个关系式之后要在 GF(2) 上解一个同样量级的稀疏线性方程组这是纯粹的计算工程问题Wiedemann 算法或者块 Lanczos 都用得上。我做的实验规模比较小用 Gaussian 消元就能过但能明显感觉到规模一上去就会卡死。LLL 则是另一条路。它的输入是一组格的基向量输出是一组更短、更接近正交的新基并且保证第一个向量不会超过最短向量的 2^((n-1)/2) 倍。这个界听起来很松但在实际应用里极为强大。Coppersmith 方法就是把求多项式的小根转成格上的短向量问题然后用 LLL 解出来从而在知道部分因子的情况下加速分解。这部分我花的时间最多因为需要转变思维方式——不再处理单个整数而是处理整数构成的向量和矩阵正交化的直觉在这里比数论定理更重要。3. 实操过程与核心环节实现3.1 环境准备与工具链选择我的实验环境最后定在 Python 3.11 gmpy2 sympy SageMath用到格和椭圆曲线时才开。这个组合的逻辑是gmpy2 负责底层大整数运算和数论原语sympy 负责符号计算和现成实现的答案校验SageMath 负责格约化和椭圆曲线这类自己写不划算的部分。fpylll 也可以用来做 LLL比 SageMath 轻量如果你不想装整个 Sage 生态可以用它。安装上没什么坑pip 装 gmpy2 时会带上 GMP 依赖Linux 下如果报编译错误装 libgmp-dev 和 libmpfr-dev 就好。SageMath 建议用官方提供的环境自己从源码编译会非常痛苦不值得。有一点必须说清楚现成库是用来对答案的不是用来交作业的。sympy 有factorint但你如果直接调它你就完全不知道 Pollard rho 的失败模式是什么遇到卡住的情况也无从下手。我的做法是自己写一版然后跟库的结果逐一比对特别是那些边界情况——平方数、2 的幂、Carmichael 数、两个大素数的乘积。3.2 Miller-Rabin 的完整实现与参数计算先看我最后定稿的版本这段代码在实验里用了很多次稳定性没问题from gmpy2 import mpz, powmod import random SMALL_PRIMES [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37] def miller_rabin(n, rounds20): n mpz(n) if n 2: return False # 小素数直接判定顺带过滤掉绝大多数偶数 for p in SMALL_PRIMES: if n % p 0: return n p # n - 1 2^s * dd 为奇数 d n - 1 s 0 while d % 2 0: d // 2 s 1 for _ in range(rounds): a random.randrange(2, int(n - 1)) x powmod(a, d, n) if x 1 or x n - 1: continue for _ in range(s - 1): x x * x % n if x n - 1: break else: return False return True几个实现上的细节值得展开。第一n - 1 2^s * d的拆解不要用位运算去数末尾零然后移位因为 Python 的整数除法在这里已经足够快而且逻辑更清楚。第二那个for...else结构是 Python 特有的else分支在循环正常结束没 break时执行正好对应整串值都没出现 n-1的情况用起来很干净但读代码的人容易懵所以我一般会加注释。第三rounds默认给 20如果你的 n 小于那个 3.3×10^24 的确定性界其实完全可以改成用固定的 13 个素数基代码会变成DETERMINISTIC_BASES [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41] def is_prime_deterministic(n): n mpz(n) if n 2: return False for p in DETERMINISTIC_BASES: if n % p 0: return n p if n 3317044064679887385961981: d, s n - 1, 0 while d % 2 0: d // 2 s 1 for a in DETERMINISTIC_BASES: x powmod(a, d, n) if x 1 or x n - 1: continue for _ in range(s - 1): x x * x % n if x n - 1: break else: return False return True return miller_rabin(n, 40)这段的意义在于在你的输入规模小于确定性界时判定结果没有概率成分实验报告里写正确是站得住的而不是高概率正确。3.3 Pollard rho 与 Brent 优化基础版 Pollard rho 我用 Floyd 判圈写代码短适合理解原理from gmpy2 import gcd, mpz def pollard_rho(n, c1): n mpz(n) if n % 2 0: return 2 x y 2 d 1 f lambda v: (v * v c) % n while d 1: x f(x) y f(f(y)) d gcd(abs(x - y), n) if d n: return None # 本轮失败换 c 重来 return d这里d n是典型的失败信号意思是撞到了 n 自己的倍数这时序列在同一圈里绕了必须换一个 c 重跑。c 的常见取值是 1、2、3……但要注意 c 0 和 c -2 会导致退化不要用。真正的性能提升来自 Brent 版本核心改动有两个一是用记录上次判圈位置的方式替代 Floyd 的双指针把每轮的计算量减半二是把 gcd 批量化。批量 gcd 的做法是维护一个累积乘积 q每 k 次迭代才做一次 gcd(q, n)k 通常取 100 左右。代码如下def pollard_rho_brent(n, c1, m128): n mpz(n) if n % 2 0: return 2 y, r, q 2, 1, 1 g 1 x ys y while g 1: x y for _ in range(r): y (y * y c) % n k 0 while k r and g 1: ys y for _ in range(min(m, r - k)): y (y * y c) % n q q * abs(x - y) % n g gcd(q, n) k m r * 2 if g n: # 回退到逐个 gcd while True: ys (ys * ys c) % n g gcd(abs(x - ys), n) if g 1: break return None if g n else g这个版本在分解 60 位左右的半素数时跟朴素 Floyd 版本相比大概有 2 到 3 倍的提速主要来自批量 gcd 和减少的取模次数。实测下来很稳。注意Pollard rho 只能一次找出一个因子要完整分解必须递归。递归的时候别忘了对拿到的因子再判一次素性否则会在素数上无限递归。3.4 BSGS 与 Pohlig-Hellman 的组装BSGS 我写成了支持指定群阶的版本因为在实际调用里阶通常是已知的有限域的阶是 p-1椭圆曲线的阶可以预先算from gmpy2 import mpz, powmod, invert def bsgs(g, h, n, orderNone): g, h, n mpz(g), mpz(h), mpz(n) N mpz(order) if order else n - 1 m mpz(int(N ** 0.5) 1) if N 10**16 else isqrt(N) 1 table {} cur mpz(1) for j in range(int(m)): table.setdefault(cur, j) cur cur * g % n factor powmod(invert(g, n), m, n) cur h for i in range(int(m) 1): if cur in table: return i * m table[cur] cur cur * factor % n return None这里table.setdefault(cur, j)是为了在小步阶段出现重复值时保留最小的 j避免解出来的 x 不是最小非负解。这个细节在验证解的时候很关键我第一次跑的时候就是因为覆盖了旧值导致返回的 x 比理论最小解大了一圈还以为是算法错了。Pohlig-Hellman 的组装逻辑是三步先把群阶 N 分解成 ∏ p_i^e_i然后对每个素数幂把问题降到 p_i 阶子群上用 BSGS 解出 x mod p_i^e_i最后用中国剩余定理把这堆同余式拼成 x。整个流程里最容易出错的是提升那一步也就是从解 x mod p 提升到 x mod p^e需要迭代地修正。我用 SageMath 的discrete_log对着写了一版校验脚本两组结果必须完全一致才算过。3.5 实验报告的写法与数据记录这门课的实验报告不要求你写漂亮话但要求数据完整。我的模板包含四块算法伪代码、关键实现、复杂度理论值与实测值对照表、失败案例记录。第四块是我觉得最有价值的因为老师看的就是你到底有没有真的跑过。下面是我做的 Pollard rho 实测表的一个片段用来感受一下增长趋势待分解 n 的位长期望迭代次数理论 n^(1/4)实测平均迭代次数平均耗时32 bit约 2^8 256310 0.01 s48 bit约 2^12 409647000.05 s64 bit约 2^16 65536780001.2 s80 bit约 2^20 ≈ 1.0×10^61.3×10^626 s实测值普遍比理论值高一点主要原因有两个一是理论界用的是渐近常数小规模下常数项影响明显二是成功找到因子需要的迭代次数是个随机变量方差很大跑十次取平均才比较可信。建议实验时每个规模至少跑 20 次记录平均值、最小值、最大值。4. 作业、实验与期末复习的组织方式4.1 作业题的三种典型类型这门课的作业基本可以归成三类识别出类型之后就知道该用什么工具。第一类是计算题比如判断以下五个数哪些是素数给出判定过程。这类题不能只写结果要写清楚你用了什么基、迭代到哪一步、为什么停下来老师是看过程给分的。第二类是证明题比如证明 Pollard p-1 中如果 p-1 是 B-光滑的则算法一定能找到因子这类题要老实按定义推。第三类是实现题给一个算法描述让你写代码并在指定数据上跑出结果。第三类最耗时间但也是最容易拿满分的因为标准明确。我的经验是把实现题做成一个可复用的小模块比如把素性检测、分解、扩展欧几里得分别写成函数放在一个numtheory.py里后面作业直接 import省下来的时间去处理那些真正难的数学题。4.2 期末复习的主线梳理期末复习我没有按章节走而是按问题—算法—复杂度这三列做了一张大表把整门课的内容压在一张纸上。表的左边是问题素性、分解、离散对数、模平方根、格上短向量中间是对应的算法右边是复杂度。这样做的最大好处是能一眼看出哪些问题是多项式时间的、哪些是亚指数的、哪些是指数的考试里但凡问可行性答案基本就在这张表上。还有几个必考的点我列一下Miller-Rabin 的误判概率推导、Pollard p-1 的光滑性条件、BSGS 的时间空间权衡、中国剩余定理在同余方程组求解中的应用、LLL 的输出向量长度上界。这几个点几乎年年出现而且都是可以在两页纸内写清楚的性价比极高。4.3 复习时容易忽略的三个细节第一个细节是位运算复杂度和算术复杂度的区别。写复杂度的时候加法是 O(log n) 而不是 O(1)乘法是 O(log²n) 而不是 O(1)因为位运算模型中数字本身要占 log n 个机器字。这一点在考试里经常被扣分我第一遍复习时也没在意后来看参考书才发现标准写法是要带 log 因子的。第二个细节是随机算法的期望复杂度和最坏复杂度。Miller-Rabin 的最坏情况可能一直抽到坏底但概率极低Pollard rho 也是期望 O(n^(1/4))最坏情况没有好界。答题时必须说清楚你说的是期望还是最坏这两者在评分里是分开的。第三个细节是素数与不可约多项式的类比。在 GF(p)[x] 里不可约多项式扮演的角色和整数里的素数完全对应所以 Miller-Rabin、Pollard rho 这些算法在多项式环里都有对应版本。课程里讲这个类比的时候我一开始没听懂后来做了一道在 GF(2)[x] 里判断不可约性的题一下就通了。5. 常见问题与排查技巧实录5.1 踩坑速查表下面这张表是我这一年攒下来的每一条都真实踩过现象大概率原因处理方式Pollard rho 卡死不返回c 取值导致序列退化或者撞上了 d n换 c 重跑加迭代上限超限就换算法Miller-Rabin 把明显合数判成素数取底用了浮点数或 n-1 拆解写错检查 randrange 的上下界打印 s 和 d 验证BSGS 内存爆掉群阶太大m √N 的表存不下改用 Pohlig-Hellman 先降阶或换 Pollard rho 求离散对数平方根判断出错用了浮点开方大整数精度丢失一律用 isqrt 或 gmpy2.irootgcd 调用耗时占比过高每轮迭代都算一次 gcd用批量 gcd累积 100 次左右再算LLL 输出向量很长格的基向量量级差太大或者参数 δ 取太小归一化输入δ 取 0.99检查格定义是否正确二次筛法解不出因子收集的关系式不够或者因子基选得太小增大因子基和光滑界 B检查线性代数部分是否真的解出了相关组合Carmichael 数通过费马测试只用了费马小定理必须升级到 Miller-Rabin费马测试不能单独用这里我想特别强调最后一条。很多人第一次写素性检测就是费马测试跑出来很快然后拿 561、1105、1729 这几个数一测就懵了。这类数在实践里一定要专门拿来做测试用例我现在的测试脚本里固定包含 {561, 1105, 1729, 2465, 2821, 6601, 8911} 这一组每次改完代码都跑一遍。5.2 性能调优与调试的几个心得关于调试我的建议是把中间量全部打出来不要用断点一行行看。计算数论的 bug 往往是数值错了而不是逻辑跳错了断点帮不上忙。我的做法是在算法里加一个debug开关打开之后打印每一轮的 x、y、gcd 结果然后对着时间序列看什么时候开始发散。Pollard rho 卡住的时候看这个日志能立刻判断是撞到 d n 了还是循环长度不对。关于性能有个反直觉的结论瓶颈经常不在你以为的地方。我一开始以为 Pollard rho 的瓶颈是模乘优化了半天模乘之后发现总耗时只降了 5%一测才发现 gcd 占了七成时间。这也是为什么批量 gcd 的收益那么大。做实验的时候一定要用cProfile先测一下热点不然就是白费劲。还有一个很实用的技巧是做规模爬坡测试。不要一上来就跑到目标规模而是从 20 位开始每次加 8 位记录耗时。这样你能观察到耗时是在哪一段开始指数上升的也能提前发现位长到某一段突然慢十倍的问题——通常是因为某处从 Karatsuba 切到了更低效的实现或者某处开始溢出到了新的算法分支。最后提一下格约化那块。LLL 的调试比前面几个算法都难因为你看不到中间数只能看到向量。我的办法是用小维度比如 3 维或 4 维手工算一遍 Gram-Schmidt把每一步的 μ 系数和新的基向量都写下来然后再跟代码的输出逐项比对。维度超过 10 之后手工就不现实了所以在能手工验证的规模上把代码调对是后面做大规模实验的前提。整理完这些笔记之后我的一个体会是USTC 计算数论这门课最有价值的地方不是让你记住多少个算法而是让你养成一种先估复杂度、再选工具、最后才动手写的习惯。这个习惯在别的领域里同样好用因为绝大多数看起来无解的问题其实是规模选错了——把规模降一个数量级或者换一个针对结构的专用算法原来跑不通的东西就会变得很轻松。我后来处理别的计算任务时脑子里第一个冒出来的问题已经变成了输入规模是多少我用的是什么复杂度这个改变就是从这门课的实验一次次超时开始的。如果你也在学这门课我的建议是别跳过任何一个看起来很简单的算法。试除法、扩展欧几里得、中国剩余定理这三个东西我在后面所有实验里都用到了写熟一遍能省下大量时间。至于那些高级的算法先把思想和复杂度搞明白实现上能跑通小规模就行真的需要大数运算的时候交给成熟的库去做把精力留给算法选择本身。