硬要说的话，佩尔方程最容易让人产生误判的地方，就是它在纸面上看着特别"初中"——x² - Dy² = 1，两个未知数，一个参数，怎么看都像是课后习题第一题。真动笔算的时候才会发现，D 取 61 这一个数，最小的正整数解就已经到了x = 1766319049、y = 226153980。我最早遇到它是在处理一个整数逼近的问题，当时想都没想就写了个循环，从 y = 1 开始往上试，试到几十万的时候觉得不对劲，一算才知道按那个速度得跑大半天。后来换成连分数，同一个答案不到一毫秒就出来了，而且全程没有出现过一个浮点数。这篇就把这套方法从头到尾拆一遍：连分数到底在做什么、√D 的展开为什么一定会循环、周期长度怎么决定你取第几个渐近分数、代码里哪几个地方最容易翻车，以及那些看着不起眼但会让程序给出错误答案的输入边界。只要你会写循环、能看懂整数除法，就能跟着跑通。
1. 那个把暴力枚举逼到墙角的方程
先把问题摆清楚。给定一个正整数 D，并且 D 不是完全平方数，求整数解x² - Dy² = 1。这个方程只要求整数解，而且一旦找到一组正整数解，就能通过代数运算生成无穷多组。所以真正的问题是：最小的一组正整数解在哪里。
为什么非要限定 D 不是完全平方数？因为如果 D = k²，方程就变成x² - k²y² = 1，也就是(x - ky)(x + ky) = 1。两个整数相乘等于 1，只有 1×1 和 (-1)×(-1) 两种可能，解出来 y = 0、x = ±1，没有正数解。所以完全平方的 D 要提前排除，这一点在写代码时必须显式判断，否则后面的连分数展开会在第一步就卡住。
1.1 它到底在问什么：一个"逼近"问题
把方程两边除以 y²，得到(x/y)² = D + 1/y²。也就是说，x/y 必须是 √D 的一个相当精确的有理逼近——它的平方只比 D 大了一点点，误差是 1/y²。
换句话说，佩尔方程和"用分数逼近无理数"是同一件事的两个说法。这一点想通了，后面的所有内容就顺了：既然要找的是 √D 的高精度有理逼近，那自然应该去找一种专门生产"最优有理逼近"的工具，而不是盲目地穷举。
这里有个常见的误解：很多人觉得 x/y 越接近 √D 越好。实际上不是"越接近越好"这么笼统，而是有一个精确的量化标准——
|x/y - √D| < 1/y²这个量级。满足这个条件的分数，才是真正有价值的逼近。
1.2 三条看起来能走、实际上走不通的路
我最初试过三种朴素思路，每一种都在 D 稍大一点的时候崩掉。
第一条：枚举 y。对每个 y 算Dy² + 1，判断是不是完全平方。D = 61 时 y 要到 226153980，也就是两亿多次循环，每次还要做一次大整数乘法和开方。即使每次只要一微秒，也是两百多秒，实际远远不止。
第二条：枚举 x。x 要到 17 亿，更慢。而且判断(x² - 1)/D是不是完全平方，同样躲不开开方运算。
第三条：用小数逼近。先算 √D 的浮点值，然后不断找分母递增的最佳逼近。这条路方向是对的，但如果不配合整数化的算法，很快就会踩到浮点精度的坑——当分母到几万、几十万的时候，双精度已经分不清相邻的两个候选分数了。
这三种思路的共同问题是：它们在"搜索"解，而佩尔方程的解不是一个能靠搜索高效定位的东西。解的位置在数轴上没有任何规律可言，D = 61 的解分母到了两亿多，D = 13 的解分母才 180，中间没有任何可预测的梯度。
2. 连分数把"无限逼近"拆成一串整数
连分数的基本形式，是把一个数写成"整数部分 + 1/(整数部分 + 1/(...))"这种嵌套。对一个实数 α，记a₀ = ⌊α⌋，然后令余项α₁ = 1/(α - a₀)，再取a₁ = ⌊α₁⌋，如此反复，就得到一个整数序列[a₀; a₁, a₂, a₃, ...]。
这个过程的妙处在于，它对有理数和无理数的行为完全不同：有理数会在有限步之后终止（因为每次取余的分数分子在严格变小），无理数则会无限延伸下去。所以连分数天然就是一个"有理数识别器"——如果某个 √D 的展开在第几步就停了，那它必然是有理数，也就意味着 D 是完全平方。
2.1 为什么连分数产生的是"最优"逼近
从连分数序列里截断，取前 n 项算成一个普通分数，得到的就叫渐近分数。这些分数有一个极其强硬的性质：在所有分母不超过 qₙ 的分数里，pₙ/qₙ 是离 α 最近的那一个。没有例外。
这个性质用生活化的说法就是：如果你的预算只能用最多 qₙ 个格子去切一条线段来逼近某个长度，那么连分数给你的切法就是最优切法，任何其他切法都不会更接近。而且它还给出一个量化的误差上界：|α - pₙ/qₙ| < 1/(qₙ · qₙ₊₁) < 1/qₙ²。
把这个不等式套回佩尔方程。设 α = √D，那么|√D - p/q| < 1/q²，两边同乘q(p + q√D)，得到：
|p² - Dq²| = |p - q√D| · |p + q√D| < (1/q) · (2q√D + 1) = 2√D + 1/q左边是个整数，右边是2√D + 1/q，也就是一个只和 D 有关、和 q 无关的常数。这就是关键：渐近分数代进p² - Dq²之后，结果的绝对值被锁死在一个上界之内，不会随着分母增大而增大。有些渐近分数会让它等于 1，有些会让它等于 -1，有些等于 2、3、4 之类的小整数，但绝不会跑到几万去。
相比之下，如果你随便拿一个分母为 q 的分数去试，p² - Dq²的量级大概是 q·√D，会跟着分母一路飙升。这就是连分数和随机枚举的本质区别。
2.2 两个递推式，全程只用整数
渐近分数的计算不需要任何浮点运算。设p₋₂ = 0, p₋₁ = 1，q₋₂ = 1, q₋₁ = 0，那么对每一项 aₙ：
pₙ = aₙ · pₙ₋₁ + pₙ₋₂ qₙ = aₙ · qₙ₋₁ + qₙ₋₂两个线性递推，全是整数加法和乘法。这一点非常重要，因为它意味着无论解涨到多少位，计算过程都不会丢精度。D = 61 的解只有十位，但有些 D 的解能到四五十位，那时候浮点早就废了，而整数递推照样精确。
我把这个递推在纸面上跑一遍 D = 2 感受一下。√2 的展开是[1; 2, 2, 2, ...]，a₀ = 1，之后全是 2。
| n | aₙ | pₙ | qₙ | pₙ² - 2qₙ² |
|---|---|---|---|---|
| 0 | 1 | 1 | 1 | -1 |
| 1 | 2 | 3 | 2 | 1 |
| 2 | 2 | 7 | 5 | -1 |
| 3 | 2 | 17 | 12 | 1 |
| 4 | 2 | 41 | 29 | -1 |
第二行就给出了 (3, 2)：9 - 8 = 1。这就是 D = 2 的最小正解。第四列那个正负交替的规律也不是巧合，后面会专门说。
顺带一提，(1, 1) 这一组满足的是
x² - 2y² = -1，不是原方程。这种差一个符号的解在连分数里出现得非常自然，处理的时候要留意别拿错了。
3. √D 的展开为什么会循环
真正让连分数能解佩尔方程的，不是"它能逼近"，而是"√D 的展开一定是循环的"，而且循环回到起点的位置可以精确判定。
3.1 三元组 (m, d, a) 的整数迭代
直接对 √D 做α → 1/(α - ⌊α⌋)这个操作，中间会引入无理数，不适合手算也不适合写代码。所以实际用的是等价的一套纯整数递推。对非完全平方的 D，令：
a₀ = ⌊√D⌋, m₀ = 0, d₀ = 1 mₖ₊₁ = dₖ · aₖ - mₖ dₖ₊₁ = (D - mₖ₊₁²) / dₖ aₖ₊₁ = ⌊(a₀ + mₖ₊₁) / dₖ₊₁⌋这三个式子里的除法都是整除，每一步的结果都是整数。这不是巧合，是可以证明的：只要初始 d₀ = 1 且 D 非完全平方，那么每一步的D - mₖ₊₁²一定能被 dₖ 整除。这个性质是整个算法能落地的基石，如果除不尽，说明输入有问题（D 是完全平方，或者中途用了浮点导致 m 算错）。
第一个式子的含义是"取余的反面"——它保证了 m 始终落在[0, √D)这个区间里。第二个式子把分母规整成新的 d。第三个式子是取整数部分，因为(a₀ + m)/d这个形式的整数部分可以直接用整数除法算出来，不需要真的去开方。
3.2 状态只有有限个，所以必然重复
m 的取值范围是0 ≤ m < √D，也就是有限多个整数。d 呢？由dₖ · dₖ₊₁ = D - mₖ₊₁²，而mₖ₊₁² < D，所以D - mₖ₊₁²是正数且小于 D，于是 d 也必然落在有限范围内。
有限个 (m, d) 组合，无限次迭代，抽屉原理直接给出结论：必然会重复。而一旦重复，后面就会原样循环下去。更进一步的结论是，第一次重复一定发生在三元组回到初始状态(m, d) = (0, 1)的时候，也就是整个循环是从 a₀ 之后严格开始的。
实际编程时不需要真的去检测"三次重复"，有一个更简洁的终止条件：当 aₖ = 2a₀ 时，循环结束。这并不是一个近似判据，而是可以严格证明的——aₖ 取到 2a₀ 的那一步，dₖ 恰好等于 1。周期长度就是那一步的下标。
我自己写第一版的时候用的是"检测 (m, d) 是否回到 (0, 1)"作为终止条件，也没问题，只是多存了一个集合。后来改成
a == 2*a0之后代码短了不少，跑了几万组测试没有出现过差异。
3.3 那个把一切串起来的恒等式
所有推导最后都汇总到一个式子上：
pₙ² - D · qₙ² = (-1)^(n+1) · dₙ₊₁这里 pₙ/qₙ 是第 n 个渐近分数，dₙ₊₁ 就是上面递推里算出来的那个 d。这个恒等式是整个方法的心脏。
它解释了三件事。第一，为什么结果被锁死——因为 d 始终落在[1, 2√D]之间，所以pₙ² - Dqₙ²的值永远在这个范围里。第二，什么时候能取到 1——当dₙ₊₁ = 1且(-1)^(n+1) = 1时，也就是 n 为奇数的时候。第三，什么时候取到 -1——同样 d = 1 但 n 为偶数。
而 d = 1 出现的时刻，正好就是循环回到起点的那些位置。所以"周期长度"这个量，直接决定了你要取第几个渐近分数才能拿到解。
4. 周期长度决定你取第几项
把上面的结论整理成可以直接用的规则。设 √D 的连分数周期为[a₁, a₂, ..., a_L]，长度是 L。
- L 为偶数：最小正解是第
L-1个渐近分数（下标从 0 开始算，a₀ 是整数部分）。 - L 为奇数：最小正解是第
2L-1个渐近分数。 - L 为奇数时，
x² - Dy² = -1有解，最小解是第L-1个渐近分数。 - L 为偶数时，
x² - Dy² = -1没有整数解。
奇偶分界这条规则的存在感非常强。为什么奇数周期只需要多跑一圈？因为每跑一圈，(-1)^(n+1)的符号就翻转一次，dₙ₊₁的序列会原样重复。所以当 L 是奇数时，第一圈末尾给出的是 -1（或者某个不等于 1 的值），第二圈末尾符号翻正，正好拿到 1。当 L 是偶数时，第一圈末尾符号已经翻正了，直接用。
举个例子体会一下周期长度和解的大小之间的关系，差别大得离谱：
| D | 连分数周期 | 周期长度 L | 最小正解 (x, y) |
|---|---|---|---|
| 2 | [2] | 1 | (3, 2) |
| 3 | [1, 2] | 2 | (2, 1) |
| 5 | [4] | 1 | (9, 4) |
| 7 | [1, 1, 1, 4] | 4 | (8, 3) |
| 13 | [1, 1, 1, 1, 6] | 5 | (649, 180) |
| 61 | [1,4,3,1,2,2,1,3,4,1,14] | 11 | (1766319049, 226153980) |
D = 7 的周期长度是 4，解只有 (8, 3)；D = 61 的周期长度是 11，解就涨到了十七亿。周期长度每增加一点，解的位数就往上跳一大截，这种增长是接近指数的。这也解释了为什么"看起来更简单的 D"反而更难算——D = 61 比 D = 62 难得多，因为 62 的周期很短。
在 D 小于一万的范围内，有些 D 的最小解 x 的位数能到四十多位。这类数用任何浮点方案都处理不了，只有纯整数递推才能给出精确结果。
5. 手推两个例子，把节奏跑顺
纸面上完整跑一遍，比看十遍公式管用。我选 D = 13 和 D = 61，一个周期短一个周期长，正好对照。
5.1 D = 13：五步跑完一个周期
a₀ = ⌊√13⌋ = 3。然后按三元组递推：
| 轮次 | m | d | a |
|---|---|---|---|
| 0 | 0 | 1 | 3 |
| 1 | 3 | 4 | 1 |
| 2 | 1 | 3 | 1 |
| 3 | 2 | 3 | 1 |
| 4 | 1 | 4 | 1 |
| 5 | 3 | 1 | 6 |
第 5 轮 a = 6 = 2a₀，同时 d = 1，循环结束。周期是[1, 1, 1, 1, 6]，长度 L = 5，奇数。
按规则，最小正解取第2L-1 = 9个渐近分数。跑递推：
| n | aₙ | pₙ | qₙ |
|---|---|---|---|
| 0 | 3 | 3 | 1 |
| 1 | 1 | 4 | 1 |
| 2 | 1 | 7 | 2 |
| 3 | 1 | 11 | 3 |
| 4 | 1 | 18 | 5 |
| 5 | 6 | 119 | 33 |
| 6 | 1 | 137 | 38 |
| 7 | 1 | 256 | 71 |
| 8 | 1 | 393 | 109 |
| 9 | 1 | 649 | 180 |
验算：649² = 421201，180² = 32400，13 × 32400 = 421200，差正好是 1。而第 4 项18² - 13 × 5² = 324 - 325 = -1，这就是负佩尔方程的最小解。两个解都在同一张表里，这也是连分数的一个便利之处。
5.2 D = 61：从十一次循环里读出十位数的解
a₀ = ⌊√61⌋ = 7。
| 轮次 | m | d | a |
|---|---|---|---|
| 0 | 0 | 1 | 7 |
| 1 | 7 | 12 | 1 |
| 2 | 5 | 3 | 4 |
| 3 | 7 | 4 | 3 |
| 4 | 5 | 9 | 1 |
| 5 | 4 | 5 | 2 |
| 6 | 6 | 5 | 2 |
| 7 | 4 | 9 | 1 |
| 8 | 5 | 4 | 3 |
| 9 | 7 | 3 | 4 |
| 10 | 5 | 12 | 1 |
| 11 | 7 | 1 | 14 |
第 11 轮 a = 14 = 2a₀，d = 1，循环结束。周期[1,4,3,1,2,2,1,3,4,1,14]，长度 L = 11，奇数。所以最小正解在第2×11-1 = 21项。
手算到第 21 项太累，我把关键节点列一下。第 10 项给出29718/3805，验算29718² - 61 × 3805² = 883159524 - 883159525 = -1，这是负方程的解。继续到第 21 项，得到1766319049/226153980，代入验算差值正好是 1。
这里能直观感受到一件事：中间那些渐近分数的pₙ² - 61qₙ²一直在 12、3、4、9、5、5、9、4、3、12 这些小数字之间来回跳，从来没有失控。这正好印证了第 3 节那个恒等式——d 的取值范围被限制在[1, 2√61]之内，所以差值不可能变大。真正长大的只有 p 和 q 本身。
6. 代码落地时真正会咬人的地方
算法本身不复杂，但要让它在各种 D 上稳定跑出正确结果，有几个地方必须处理干净。我按自己踩坑的顺序说。
6.1 展开函数：把状态压缩成三个整数
from math import isqrt def cf_sqrt(D): """返回 (a0, period)，period 为 sqrt(D) 的连分数周期列表。 D 为完全平方时返回 (a0, [])。""" 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 a = (a0 + m) // d period.append(a) if a == 2 * a0: return a0, period这段代码里最值得注意的是isqrt。它是整数平方根，返回的是精确的 ⌊√D⌋，不会有任何舍入误差。int(math.sqrt(D))在 D 接近 10¹⁸ 的时候就会开始出错，因为 float 的尾数只有 53 位，超过 2⁵³ 之后相邻整数已经分不开了。
还有一个细节是//整除。这里的除法在数学上是精确整除，用//和/结果一样，但//会保持整数类型，避免意外产生浮点。如果哪天除不尽了，//会悄悄丢掉余数而不是报错，所以调试阶段可以加一行断言：
assert (D - m * m) % d == 0我加过这行，跑了几千个 D 都没触发过。一旦触发，基本可以确定是 D 是完全平方却漏判了。
6.2 取第几项：把周期拼成足够长的序列
拿到 (a₀, period) 之后，需要先确定目标下标 n，再把 a 序列补齐到 n 位。
def pell_min(D): """返回 x^2 - D*y^2 = 1 的最小正整数解 (x, y)。 若 D 为完全平方则返回 None。""" a0, period = cf_sqrt(D) if not period: return None L = len(period) n = (L - 1) if L % 2 == 0 else (2 * L - 1) seq = [a0] while len(seq) <= n: seq.extend(period) p_prev2, p_prev1 = 0, 1 q_prev2, q_prev1 = 1, 0 for i in range(n + 1): a = seq[i] p = a * p_prev1 + p_prev2 q = a * q_prev1 + q_prev2 p_prev2, p_prev1 = p_prev1, p q_prev2, q_prev1 = q_prev1, q return p_prev1, q_prev1跑一下几个已知值验证：
for D in [2, 3, 5, 7, 13, 61, 109]: print(D, pell_min(D))输出应该是2 (3, 2)、3 (2, 1)、5 (9, 4)、7 (8, 3)、13 (649, 180)、61 (1766319049, 226153980)、109 (158070671986249, 15140424455100)。D = 109 那组解已经到十五位了，正好用来确认大整数路径没出问题。
补 a 序列的时候我一开始写的是
seq = [a0] + period * 3，对小 D 凑合能用，但周期长度超过 3 的时候会越界。改成 while 循环补齐就没有这个问题，多跑不了几次循环。
6.3 三个容易翻车的实操细节
第一，判完全平方别用浮点。这是最经典的一个坑。写r = int(math.sqrt(n)); return r*r == n看起来很干净，但 n 一旦超过 2⁵³，math.sqrt的返回值可能偏高，int()截断之后 r 就大了 1，判断直接失效。正确的是r = isqrt(n); return r*r == n。这个坑在 D 比较小的时候完全看不出来，等测试数据涨到十几位才突然冒出来，非常难查。
第二，别在大数上做乘方再取模。有些实现会用x² mod D来判断，理论上可行，但中间结果可能超出机器字长。语言层面如果没有自动大整数（比如 C 或 Go 的固定宽度整型），到 D = 61 这一级别就已经溢出了。要么用大整数库，要么确保每一步都在 64 位以内——但 D = 61 的解需要 61 位左右的中间量，勉强卡在边界上，再大一点的 D 就必然溢出。
第三，注意下标从 0 开始还是从 1 开始。这个坑纯粹是自找的，但确实容易犯。恒等式pₙ² - Dqₙ² = (-1)^(n+1) · dₙ₊₁里，n 是从 0 开始数的渐近分数序号，而dₙ₊₁用的是递推下一轮的 d。如果代码里把这两个下标搞混，结果会在"对"和"差一个符号"之间反复横跳，而且小 D 上可能碰巧还是对的，特别迷惑。
7. 负佩尔方程与那些必须提前排除的输入
x² - Dy² = -1这个方程经常和原方程一起出现，因为连分数顺手就能给出它的解。但它不是永远有解，判断依据就一条：周期长度 L 为奇数时有解，为偶数时无解。
道理和前面一样，(-1)^(n+1)在每一圈末尾翻转一次符号，d 则原样重复。奇数周期意味着第一圈末尾的符号是 -1，此时d = 1就给出p² - Dq² = -1。偶数周期第一圈末尾符号是 +1，永远碰不到 -1 那个分支。
def pell_neg(D): """返回 x^2 - D*y^2 = -1 的最小正解，无解返回 None。""" a0, period = cf_sqrt(D) if not period: return None L = len(period) if L % 2 == 0: return None n = L - 1 seq = [a0] while len(seq) <= n: seq.extend(period) p_prev2, p_prev1 = 0, 1 q_prev2, q_prev1 = 1, 0 for i in range(n + 1): a = seq[i] p = a * p_prev1 + p_prev2 q = a * q_prev1 + q_prev2 p_prev2, p_prev1 = p_prev1, p q_prev2, q_prev1 = q_prev1, q return p_prev1, q_prev1D = 13 时返回 (18, 5)，D = 61 时返回 (29718, 3805)，D = 3 时返回 None。三个结果都对得上手算。
7.1 三种必须显式挡掉的输入
D 是完全平方。前面说过，这种情况下原方程只有 y = 0 的平凡解，负方程更是完全没有解。连分数展开在第一步d₀ = 1之后就会算出d₁ = D - a₀² = 0，再往下除零直接崩溃。所以a0 * a0 == D这个检查必须放在最前面。
D = 1。严格来说 1 是完全平方，会被上一条挡掉。但如果代码里图省事用了D > 1就往下走，D = 1 会让递推逻辑错乱。这个属于边界值处理，顺手加个判断就行。
D 含平方因子。比如 D = 12 = 4×3。这时候方程x² - 12y² = 1有解，但最小解没法直接套"周期长度"规则，因为 √12 = 2√3 的连分数展开和有理性混在一起了。正确的处理是先把平方因子提出来，转化成x² - 3(2y)² = 1，也就是解u² - 3v² = 1且要求 v 是偶数，解出来再换回 x = u、y = v/2。如果 v 不是偶数，就取下一组解（因为解是无穷多的，总能找到 v 为偶数的）。
这一条我在实际使用中吃过亏。当时输入是一批随机生成的 D，里面混着 8、12、18 这种含平方因子的数，程序没报错，返回的结果代入验算却差了 4 倍。查了半天才发现是没做平方因子剥离。现在我的做法是包一层预处理，先把 D 化成
D' · s²的形式，D' 无平方因子，再把解映射回去。
7.2 无穷多解怎么生成
找到最小解(x₁, y₁)之后，剩下的解可以用代数方式一次全生成出来。核心观察是：如果(x₁, y₁)是解，那么(x₁ + y₁√D)ⁿ展开之后的有理部分和 √D 系数部分，也一定构成解。
具体递推式是：
xₙ₊₁ = x₁ · xₙ + D · y₁ · yₙ yₙ₊₁ = x₁ · yₙ + y₁ · xₙ拿 D = 2 验证，最小解是 (3, 2)。下一组：x₂ = 3×3 + 2×2×2 = 17，y₂ = 3×2 + 2×3 = 12。验算17² - 2×12² = 289 - 288 = 1。再下一组 x₃ = 3×17 + 2×2×12 = 99，y₃ = 3×12 + 2×17 = 70，`99² - 2×70² = 980
