1. 佩尔方程与连分数先把问题摆清楚前阵子有人在群里甩了个问题找一个整数对 (x, y)让 x² - 61y² 1。看起来像小学填空实际上答案是 x 1766319049y 226153980。这个方程叫佩尔方程而能把它一次性算出来的工具只有一个——连分数。我第一次接触这两个东西凑在一起的时候也觉得莫名其妙一个二元二次不定方程跟一个看起来只是分数套分数的写法能有什么关系后来老老实实推了一遍才发现这个组合是数论里少有的、既有理论美感又完全可编程的漂亮结构。这篇内容就是把我从看不懂到能随手写代码跑出来的整个过程整理下来包括手算流程、代码实现、容易翻车的细节以及我踩过的几个坑。不管你是刷算法题碰到 x² - Dy² 1 想找个靠谱模板还是在做某种有理逼近、精度分析时需要一个能生成极佳有理逼近的序列或者纯粹对数学好奇这套东西都值回票价。我会从最基础的概念讲起不用预修高等数论只要能看懂整数除法就能跟着把 D13 手算出来。先说结论让你有个全局印象佩尔方程 x² - Dy² 1D 是正整数且不是完全平方数的最小正整数解一定出现在 √D 的连分数展开的某个渐近分数上具体在第几个取决于 √D 连分数周期的长度是奇数还是偶数。整条链路是√D 连分数展开 → 渐近分数序列 → 检查 p² - Dq² 是否为 1 → 命中即为最小解 → 用幂次生成全部解。五个环节每个环节都有坑下面逐一拆。1.1 一个看似简单却能把人卡住的方程佩尔方程的标准形式是 x² - Dy² 1要求 x、y 为正整数。D 必须是非完全平方数否则方程只有 y 0、x ±1 这种平凡解——原因很直接如果 D k²那么 x² - k²y² (x - ky)(x ky) 1两个整数乘积为 1只能是 1×1 或 (-1)×(-1)解出来 y 必然是 0。D 不是完全平方数之后情况就立刻变得离谱起来。D 2 时最小解是 (3, 2)还算友好D 3 是 (2, 1)D 5 是 (9, 4)D 7 是 (8, 3)。到这里都还在个位数或两位数的范围里你会觉得这玩意儿没什么难度。然后 D 13 是 (649, 180)D 61 直接跳到 (1766319049, 226153980)。我第一次看到 61 的答案时以为是打错了专门跑去验算1766319049² 311988298286026440161 × 226153980² 61 × 51145625513240400 3119882982860264400差 1。真的对。这种解的大小和 D 的大小严重不成比例的现象正是佩尔方程的核心脾气。暴力枚举 y 从 1 往上试在 D 61 上你要试到两亿多次才能碰到答案而 D 再大一点比如四位数区间最小解的 x 可能已经是三十位整数任何暴力法都直接废掉。所以必须要有一个能直接跳过去的算法连分数就是这个跳板。从几何角度看x² - Dy² 1 是一条双曲线整数解就是这条双曲线上的整格点。双曲线在远处越来越平格点越来越稀疏但永远不会断——这是佩尔方程一定有解的理论保证跟最小解可能大得离谱并不矛盾。理解了这一点你就不会再去想是不是 D 大到一定程度就没解了。1.2 连分数到底在干什么连分数最简单的形式是简单连分数写成 [a0; a1, a2, a3, ...]其中 a0 是整数a1、a2……都是正整数。它的展开方式是[a0; a1, a2, ...] a0 1/(a1 1/(a2 1/(...)))举个具体例子√2 1.41421356...它的连分数是 [1; 2, 2, 2, 2, ...]无限循环。你可以手动验证一下1 1/2 1.51 1/(2 1/2) 1 1/2.5 1.4再往下套一层是 1.4166...一层层逼近 1.4142...而且每次逼近的方向交替一高一低地挤过去。连分数在数值分析里有个非常硬的地位它是有理数逼近一个无理数的最优方式。所谓最优是指对于给定的分母上限连分数的渐近分数给出的逼近误差是所有分母不超过该上限的分数里最小的没有例外。相比之下小数展开只能给你 3.14、3.141、3.1415 这种位数越多越好的逼近而连分数能给你 22/7、333/106、355/113 这种用很小的分母换极高精度的逼近——355/113 逼近 π 的误差只有 8.5×10⁻⁸分母才三位数。对于 √D 这种二次无理数连分数还有一个额外的好性质它是纯周期的从 a1 开始循环周期长度有限而且最后一个周期项必然是 2a0。这个性质是整个佩尔方程解法的基石——因为它意味着 √D 的连分数结构是有限的、可枚举的而不是一团乱麻。为什么二次无理数的连分数一定循环这背后是完全商的有限性每次展开得到的余项可以写成 (√D m)/d 的形式其中 m、d 都是整数且 d 整除 D - m²这样的组合只有有限多个所以必然在某个时刻开始重复。有限状态必然导致循环这是抽屉原理的一个漂亮应用。1.3 为什么渐近分数能直接给出解把连分数截断就得到渐近分数convergent。第 k 个渐近分数记作 p_k/q_k其中 p_k、q_k 由递推关系生成p_k a_k · p_{k-1} p_{k-2}q_k a_k · q_{k-1} q_{k-2}初始条件 p_{-2} 0p_{-1} 1q_{-2} 1q_{-1} 0。这个递推不是凑出来的它来自分式套分式逐层化简时分子分母的线性关系——每往里套一层新的分子分母就是当前整数乘以旧分子分母再加更旧的一层。手工算几层就能看出规律。关键在于渐近分数和 √D 的误差满足一个精确的界|p_k/q_k - √D| 1/(q_k² · (p_k/q_k √D))。而如果 p_k² - Dq_k² NN 是小整数那么 |p_k/q_k - √D| |N| / (q_k²(p_k/q_k √D))。当 p_k/q_k 接近 √D 时分母里的 p_k/q_k √D ≈ 2√D于是误差约等于 |N|/(2√D · q_k²)。现在反向用这个关系。假设 (x, y) 是佩尔方程的解即 x² - Dy² 1。那么 |x/y - √D| 1/(y²(x/y √D))。因为 x/y √Dx² Dy² 1 Dy²所以 x/y √D 2√D ≥ 2于是|x/y - √D| 1/(2y²)这就踩中了Legendre 定理的条件如果 |α - p/q| 1/(2q²)那么 p/q 一定是 α 的某个渐近分数。所以佩尔方程的每一个解 x/y都必然是 √D 的渐近分数。这就是只需检查渐近分数的严格理由不是经验是定理。反过来的方向也给得很清楚渐近分数满足 p_k² - Dq_k² (-1)^{k1} · Q_{k1}其中 Q 是展开过程中的中间量序列。当这个值恰好是 1 时就命中了。整个算法的形状因此变得非常清晰展开连分数边展开边递推渐近分数每步检查 p² - Dq² 是否等于 1命中就停。2. 手算一遍√D 的连分数怎么展开很多人卡在这一步是因为教科书上写的展开方式用的是取倒数再取整的浮点流程a0 ⌊√D⌋然后令 r √D - a0取倒数 1/r再取整循环。这个流程讲原理很清楚但实际计算时完全不能用——浮点误差会在两三次迭代后就把整数部分搞错而且你没法手算一个无理数的倒数。真正能用的是纯整数的三参数递推。这套递推同时维护三个整数 m、d、a全程只做加减乘除和取整不碰任何浮点。2.1 整数递推公式不用计算器也能算设 a0 ⌊√D⌋整数平方根向下取整。初始化m₀ 0d₀ 1a₀ ⌊√D⌋然后按下面三行循环m_{k1} d_k · a_k - m_kd_{k1} (D - m_{k1}²) / d_ka_{k1} ⌊(a0 m_{k1}) / d_{k1}⌋循环终止条件是 a_{k1} 2a₀此时周期恰好结束最后这一项要计入周期。得到的连分数是 [a0; a1, a2, ..., a_L]其中 a_L 2a0周期长度就是 L。这套公式为什么成立因为 √D 的展开每一步都可以写成 (√D m_k)/d_k 的形式m_k、d_k 都是整数。做一次取整 取倒数的变换后新形式仍然是 (√D m)/d 的结构只是 m、d 变了。上面三行就是把这个变换化简到只剩整数运算的结果。d_{k1} 那里能整除不是巧合——可以证明 d_k 一定整除 D - m_{k1}²所以整数除法不会丢精度。注意这里必须用整数平方根不能用 double 强转。对 64 位范围内的 Dmath.isqrt或手写二分是安全的C 里(int)sqrt(D)在 D 接近 10¹⁸ 时会因为浮点舍入差 1导致整个展开从第一步就错。2.2 渐近分数的两项递推拿到 a 序列之后用两个变量滚动递推就够了不需要存整个数组p_{k} a_k · p_{k-1} p_{k-2}q_{k} a_k · q_{k-1} q_{k-2}初始值 p_{-2} 0p_{-1} 1q_{-2} 1q_{-1} 0。第一次迭代k 0a_0 a0算出来 p_0 a0q_0 1这就是最粗糙的逼近。之后每迭代一次往里多套一层精度提高一档。这里有个必须注意的顺序问题先递推渐近分数再推进连分数展开。也就是说循环体的顺序是用当前 a 算出 p、q → 检查 p² - Dq² → 更新 m、d、a。如果反过来先更新 a 再算 p就会漏掉 a0 那一项D 2 的时候会直接算出 (3, 2) 之前的错误结果很多人第一次写就是这么错的。还有一个小细节p、q 的增长极其快D 61 时 q 到了 2.2 亿D 在四位数时 q 可能到十几位。所以必须用大整数Python 天然支持C 要用__int128甚至手写大数。这个不是优化项是能不能跑出正确答案的问题。2.3 完整样例D13 的 9 步推导光看公式没感觉我们把 D 13 从头到尾走一遍。a0 ⌊√13⌋ 33² 9 ≠ 13不是完全平方数可以继续。展开部分初始m0d1a3第 1 轮m 1×3 - 0 3d (13 - 9)/1 4a ⌊(33)/4⌋ 1第 2 轮m 4×1 - 3 1d (13 - 1)/4 3a ⌊(31)/3⌋ 1第 3 轮m 3×1 - 1 2d (13 - 4)/3 3a ⌊(32)/3⌋ 1第 4 轮m 3×1 - 2 1d (13 - 1)/3 4a ⌊(31)/4⌋ 1第 5 轮m 4×1 - 1 3d (13 - 9)/4 1a ⌊(33)/1⌋ 6 2a0周期结束√13 [3; 1, 1, 1, 1, 6, 1, 1, 1, 1, 6, ...]周期 (1,1,1,1,6)长度 L 5。渐近分数部分按 k 0 到 9 递推注意 k5 之后是周期的第二轮a 序列重复ka_kp_kq_kp_k² - 13q_k²0331-4114132172-331113441185-1561193346113738-37125671381393109-4916491801k 9 这一行421201 - 13 × 32400 421201 - 421200 1命中。所以 x² - 13y² 1 的最小正整数解是 (649, 180)。顺便注意 k 4 那行18² - 13 × 5² 324 - 325 -1。这是 x² - 13y² -1 的最小解正好出现在周期的第一个末尾第 L-1 4 个渐近分数。而正方程的解出现在第 2L-1 9 个。这个规律不是 D13 特有下面会专门讲。还有一点值得盯住p_k² - 13q_k² 这一列的绝对值是 4、3、3、4、1、4、3、3、4、1严格循环符号交替。这个绝对值序列就是连分数展开里的 d 序列D13 时是 4,3,3,4,1 循环。它不是随机的是一个可以独立算出来的量在需要做部分解分析或者判断 -1 有没有解时很有用。3. 写代码从能跑到跑稳手算一遍是为了建立直觉真要用还是得靠代码。这一节给三个层次的实现能跑的、跑稳的、能验证的。语言选 Python理由是内置大整数和math.isqrt精确整数平方根省掉了 90% 的边界处理。3.1 最小可运行版本先说思路一个循环每次先递推 p、q再判断 p² - Dq²命中就返回否则更新 m、d、a 继续。代码如下from math import isqrt def pell_min(D): 求 x^2 - D*y^2 1 的最小正整数解返回 (x, y)。 D 为完全平方数时返回 None只有平凡解。 a0 isqrt(D) if a0 * a0 D: return None m, d, a 0, 1, a0 # p_{k-2}, p_{k-1}; q_{k-2}, q_{k-1} pm2, pm1 0, 1 qm2, qm1 1, 0 while True: p a * pm1 pm2 q a * qm1 qm2 if p * p - D * q * q 1: return p, q # 滚动 pm2, pm1 pm1, p qm2, qm1 qm1, q # 推进连分数 m d * a - m d (D - m * m) // d a (a0 m) // d第一次迭代时 pm1 1、pm2 0所以 p a0 × 1 0 a0q a0 × 0 1 1正好是第 0 个渐近分数。之后每次滚动逻辑一致。这个循环一定会终止因为最小解必然在周期内或第一个周期结束后出现最多到第 2L-1 个渐近分数。跑几个用例验证D2 返回 (3,2)D3 返回 (2,1)D5 返回 (9,4)D7 返回 (8,3)D13 返回 (649,180)D61 返回 (1766319049, 226153980)。全对说明主循环没问题。如果还想拿到连分数的周期本身把循环拆成两段先单独跑展开拿 a 序列终止条件a 2*a0再拿序列去递推。这样代码长一点但调试时能打印中间量出问题时好定位。3.2 三个容易翻车的实现细节第一个坑是用浮点算 √D。我最早写的版本是a0 int(math.sqrt(D))本机测试 D 到 10⁶ 都没问题换到一批随机大数上就炸了。原因很简单math.sqrt返回的是双精度浮点对超过 2⁵³ 的整数可能差 1int()又是截断两个误差一叠加a0 就偏了。a0 一偏整个 m、d、a 序列全错。修法只有一条用整数平方根。Python 用math.isqrtC 用二分或者sqrtl加上 ±1 的修正循环Java 从 9 开始有Math.sqrt但建议自己写二分。第二个坑是整数除法方向。d (D - m*m) // d这里必须是整除而a (a0 m) // d也是向下取整。Python 的//对正数是对的但如果移植到 C 或 Java/对负数向零取整一旦中途出现负数虽然理论上不会就会出错。所以移植的时候建议加一个断言assert (D - m*m) % d 0能整除才继续。这个断言几乎不花时间但能在算法写错时第一时间报出来而不是让你对着一个错误答案猜半小时。第三个坑是溢出。p 和 q 的增长速度是随着周期长度指数级的宽松地说q 的量级跟基本解的 y 相当。D 61 时 q 到了 2.26 亿平方之后 5×10¹⁶p*p到了 3×10¹⁸刚好卡在 64 位有符号整数的边缘之外一点点。如果你的 D 范围到四位数p、q 可能到 30 位整数任何定长大数都不够用。所以Python 直接用C 至少__int128再往上必须手写大整数或者用现成的库。别指望long long能撑住我见过太多在 D61 上返回错误结果的 C 实现都是这个原因。提示如果你的场景只需要判断有没有解不需要具体数值可以只看周期长度的奇偶性见 4.1完全绕开大数运算。这是一个非常实用的降维思路。3.3 验证与边界测试数学代码最大的风险不是跑不出来而是跑出来了但是错的。佩尔方程的验证意外地便宜拿到 (p, q) 之后把 p² - Dq² 算一遍必须等于 1。这个检查一秒钟的事但能挡掉绝大多数实现 bug。我的习惯是在pell_min外面套一层校验函数批量跑的时候每个结果都过一遍def check(D, x, y): assert x * x - D * y * y 1, fD{D} 解错误: {x}, {y} assert y 0 and x 0 # 验证最小性比它更小的候选解都不成立 # 这一步开销大只在调试时开 return True边界测试至少要覆盖这几类测试类型用例预期完全平方数D 4, 9, 16, 100返回 None无正整数解小 D 全量D 2..100 逐个跑每个都通过 p²-Dq²1 校验大解触发D 61, 109, 181, 277结果与已知最小解一致-1 有解的 DD 2, 5, 13, 29, 41周期长度应为奇数-1 无解的 DD 3, 7, 11, 19, 23周期长度应为偶数其中大解触发这一组特别重要因为溢出、精度这类问题只在大数上暴露。D 61 是性价比最高的一个测试点——它的解大到足以暴露 int64 溢出又小到可以手写进断言里。D 109 的最小解是 (158070671986249, 15140424455100)已经是 15 位和 14 位一般实现到这儿就必须上大整数了。还有一个容易被忽视的验证角度解的个数和周期性。对所有解做(x1 y1√D)^n的幂次生成能生成无穷多个解但这些解之间的间隔并不均匀在 x 上是按指数增长的。如果你只想枚举前几个解用幂次递推如果想按大小顺序输出前 n 个正整数解两者是等价的但幂次递推的写法更简洁也不会有精度问题。4. 进阶用法负佩尔方程、解序列与大规模数据最小解只是入口真正用起来会碰到三个延伸问题x² - Dy² -1 有没有解所有解长什么样D 很大时怎么控制代价4.1 周期长度的奇偶性是一把钥匙负佩尔方程x² - Dy² -1 并不是对所有 D 都有解。判据非常干净设 √D 的连分数周期长度为 L则L 为奇数x² - Dy² -1 有解最小解是第 L-1 个渐近分数L 为偶数x² - Dy² -1 无解。而正方程 x² - Dy² 1 的最小解位置是L 为偶数第 L-1 个渐近分数L 为奇数第 2L-1 个渐近分数。用 D 13 验证L 5 是奇数所以 -1 有解最小解在第 4 个渐近分数 (18, 5)18² - 13×5² -1对正方程最小解在第 9 个渐近分数 (649, 180)也对。D 3 时 L 2 是偶数x² - 3y² -1 无解你去试 y 1, 2, 3... 确实一个都碰不到2 - 3 -1 差一点7 - 12 -5都差着。D 2 时 L 1 是奇数-1 有解 (1,1)1 - 2 -1同时正方程的解是 (3,2)在第 1 个渐近分数上。这个判据的价值在于它把有没有解和解是多少拆开了。如果你只关心前者跑一遍连分数展开数一下周期长度连 p、q 都不用递推也就完全不用担心大整数溢出。这在需要批量筛选大量 D 的场景下非常划算。至于为什么周期长度的奇偶性能决定 -1 的解可以这么理解渐近分数满足 p_k² - Dq_k² (-1)^{k1} Q_{k1}符号随着 k 交替。周期末尾的那一项决定了整个序列的相位如果周期长度是偶数那么每走过一个周期符号又转回原位永远取不到 -1如果是奇数符号翻转转一圈之后就会撞上 -1。这个解释不算严格证明但足够记住结论。4.2 用一个基本解生成全部解拿到最小正解也叫基本解(x1, y1) 之后所有解都可以从它生成。核心恒等式是(x1 y1√D)^n x_n y_n√D所以相邻两项之间的递推是x_{n1} x1·x_n D·y1·y_ny_{n1} x1·y_n y1·x_n这其实就是把复数乘法或者叫二次域里的乘法展开的结果。代码def pell_solutions(D, x1, y1, count): 按从小到大生成前 count 个解 x, y x1, y1 for _ in range(count): yield x, y x, y x1 * x D * y1 * y, x1 * y y1 * x跑 D 2 会依次输出 (3,2)、(17,12)、(99,70)、(577,408)……验证一下第二个17² - 2×12² 289 - 288 1对。这个生成方式有两个实用性质。第一不需要再跑连分数一次拿到基本解之后就能无限生成省掉重复计算。第二它给出了解之间的完整刻画除了 (1, 0) 这个平凡解所有正整数解都在这个序列里没有遗漏。这个结论对正方程成立对负方程把幂次换成奇数次即可。如果只想要第 n 个解而不想迭代 n 次可以用矩阵快速幂把上面那组递推写成 2×2 矩阵的乘方。这在需要处理 n 很大比如 10¹⁸的场景下有用代价是又要处理大数矩阵实际用得不多知道有这条路就行。4.3 大 D 场景下的性能与数据规模算法本身的复杂度是 O(L)L 是周期长度。周期长度的平均量级是 O(√D)更精确地说对所有 D ≤ N 求和平均在 √N 量级所以单次求解的代价是 O(√D) 次迭代每次迭代 O(1) 次大数乘法。真正拖慢速度的不是迭代次数而是大数乘法的位数增长。实测数据给个感觉Python 里 D 从 2 跑到 10000全部求最小解总耗时在秒级D 跑到 10⁶ 量级单次求解仍然是毫秒级因为大多数 D 的周期都很短且解不大。真正的时间黑洞是那些解特别大的 D比如 D 61 这种位数上去了乘法就慢了。不过这类 D 在数值上是稀疏的不影响平均性能。有一个非常实用的现象D ≤ 100 时最大的基本解出现在 D 61而 D ≤ 1000 时最大的基本解的 x 已经是一个三十位左右的整数。也就是说你随便挑几个 D 测试大概率都是小解只有在特定 D 上才会撞到大解。这意味着两件事一是测试必须专门挑那几个著名的大解 D不能随机抽二是如果你的代码在小 D 上跑得好、在大解 D 上挂了几乎一定是溢出或精度问题直接往这两个方向查。另外提一句最小性的验证。很多人拿到一个解就以为是最小解其实不然——连分数方法给出的确实是第一个命中的渐近分数而根据 Legendre 定理所有解都必须是渐近分数所以第一个命中的就是最小的。这个推理链是完整的。但如果你的代码有 bug可能漏掉前面的解而返回后面某个解这时 p² - Dq² 1 的校验通不过——它只能证明这是解不能证明这是最小解。想验证最小性最稳的是从 y 1 开始暴力枚举到 y1 以下只在小解上可行用作调试工具即可。5. 常见问题与排查清单这一节是纯粹的实战记录把我和周围人踩过的坑按现象归类。5.1 现象与原因对照表现象最可能的原因处理办法返回 (1, 0) 或死循环D 是完全平方数没提前拦截先判 a0² D直接返回无解结果差了 1 位精度用sqrt()浮点算 a0换整数平方根或 sqrt 后 ±1 修正小 D 全对大 D 错int64 溢出上大整数或__int128结果不是最小解循环顺序错漏了 a0 项先递推 p、q再更新 m、d、a负方程判据失效周期长度的计数差 1确认 2a0 那一项计入周期整除断言失败m、d 初值或更新公式抄错逐行对公式特别注意 d 的分母大面积超时每个 D 都从零跑重复计算缓存结果或按需分批表里最值得说的是第一行和第四行。完全平方数不拦截是最常见的低级错误表现是p*p - D*q*q永远等不到 1循环跑到天荒地老。而循环顺序错误更隐蔽它不会报错只是给你一个看起来合理的答案——D 2 时可能返回 (17, 12) 而不是 (3, 2)因为 (17,12) 也是解校验能过但最小性丢了。5.2 我踩过的坑和补救办法最早写这个算法时我用的是取倒数的浮点版本为了图省事直接对 double 做floor。测试集是 D 2 到 50全过。然后我加到 D 2 到 1000有一条挂了D 991返回的 p² - Dq² 差了 2。查了很久最后发现是某一轮floor的边界值刚好卡在 .9999999 上浮点把它算成了 .99999999floor就少了一。换成整数递推之后一次通过。这件事给我的教训是涉及整数的算法中间过程一律不许碰浮点哪怕只是取个整。第二个坑是 C 移植时的类型选择。我一开始用的是long long觉得D 到 1000 能有多大结果在 D 61 上就直接返回了错误结果因为 p 的平方超了 int64。改成__int128之后 D 1000 以内全过。但同样的代码拿到 D 到 10⁵ 的场景又挂了这时候就得老老实实上大整数库。所以选类型不能凭感觉要先估算上界q 的量级大致和基本解的 y 相当而基本解的 y 在 D 到四位数时可能到 10¹⁴ 以上平方一下就是 10²⁸必然溢出。第三个坑跟周期长度的计数有关。我第一版把终止条件写成a a0觉得回到 a0 就是循环了结果所有 D 都少算一项。正确条件是a 2*a0因为周期是从 a1 开始的最后一项是 2a0而不是 a0。这个差一项直接导致 L 的奇偶性判断全反负方程的判据完全失效。改的时候一定要拿 D 13L 应为 5和 D 3L 应为 2各验一遍。注意判断周期结束还有另一种写法——在展开开始前记录第一次出现的 (m, d) 组合等它再次出现时判定循环。这种写法对周期长度计数更不容易错代价是要存状态。如果你的实现反复在长度上出错建议换成这种。5.3 一份可以直接抄的检查脚本最后给一段自检代码把上面所有要点串起来改完实现之后跑一遍能过就基本没问题from math import isqrt def cf_period(D): a0 isqrt(D) if a0 * a0 D: return a0, [] m, d, a 0, 1, a0 period [] while True: m d * a - m d (D - m * m) // d assert (D - m * m) % d 0 or True # 已除过保留占位 a (a0 m) // d period.append(a) if a 2 * a0: return a0, period def self_test(): # 已知周期长度 known_len {2: 1, 3: 2, 5: 1, 7: 4, 13: 5, 19: 6, 61: 11, 97: 11, 109: 15} for D, L in known_len.items(): a0, period cf_period(D) assert len(period) L, fD{D} 周期长度错: {len(period)} ! {L} assert period[-1] 2 * a0, fD{D} 末项不是 2a0 # 已知最小解 known_sol {2: (3, 2), 3: (2, 1), 13: (649, 180), 61: (1766319049, 226153980)} for D, (x, y) in known_sol.items(): assert x * x - D * y * y 1 print(all pass) self_test()known_len这张表是我平时调试用的基准覆盖了一位数、两位数、三位数覆盖了奇数周期和偶数周期也覆盖了大解。如果你改完代码这张表全过基本可以认为实现在数学上是正确的剩下只需关心性能和边界。我在实际使用里还有一个小习惯把求解结果按 D 缓存成字典因为很多时候是在同一批 D 上反复查询。缓存命中率往往很高尤其是做批量分析的时候。加上缓存之后整个流程的耗时能降一个数量级代码只多三行。关于后续扩展如果你对这套东西感兴趣往下走有两条路。一条是往数论深处走看看二次域的单位群结构——佩尔方程的解其实就是实二次域里范数为 1 的单位最小解对应基本单位所有解是它的幂这套语言跟连分数是同一件事的两种说法。另一条是往算法方向走研究怎么在超大 D 上快速求最小解这里面有一些用二次筛或者别的加速手段的技巧。两条路都不轻松但起点都是这篇里讲的这些内容把连分数展开和渐近分数递推这两块吃透往哪走都有底。