类Pell方程的一般解法

2026/04/06

在Project Euler中，连续遇到了很多问题涉及到正整数变量的二元二次方程，它们一般可以表示为类Pell方程： \[x^2 - Dy^2 = N.\] 其中\(D>0\)。（当\(N=1\)时，为标准Pell方程。）我请教了Gemini，它帮我总结了这类方程的一般解法。

\(D\) 为完全平方数的情况

如果\(D=k^2\)，则方程转化为 \[(x-ky)(x+ky)=N.\]

解法为：

  1. 枚举\(N\)的所有正负因子对\((d_1,d_2)\)，使得\(d_1d_2 = N\)。
  2. 解方程组\(\begin{cases}x-ky=d_1,\\x+ky=d_2.\end{cases}\)
  3. 筛选出正整数解。

如果\(D\)不是完全平方数，则需要用下面的求解步骤。

求解标准Pell方程\(u^2 - Dv^2 = 1\)

当\(D\)不是完全平方数时，首先需找到\(u^2 - Dv^2 = 1\)的最小正整数解\((u_1, v_1)\)。

步骤：

  1. 将\(\sqrt D\)展开为连分数形式：\(\sqrt D= [a_0; \overline{a_1, a_2, \dots, a_n, 2a_0}]\)。
  2. 计算该连分数的渐近分数\(\frac{p_k}{q_k}\)。
  3. 第一个满足\(p_k^2 - D q_k^2 = 1\)的\((p_k, q_k)\)即为\((u_1, v_1)\)。

寻找类Pell方程的基础解

方程\(x^2 - Dy^2 = N\)的解分布在若干个“解类”中。每一类由一个基础解\((x_i, y_i)\)代表。

方法：枚举\(y\in[0,y_\max]\)，其中上界 \[y_\max = \begin{cases} \frac{v_1\sqrt N}{\sqrt{2(u_1+1)}}&N>0,\\ \frac{v_1\sqrt{|N|}}{\sqrt{2(u_1-1)}}&N<0. \end{cases}\] 对每个\(y\)，如果\(Dy^2 + N\)是完全平方数，则\(x=\sqrt{Dy^2 + N}\)是整数。此时\((x, y)\)是一组基础解。

构造通解

一旦找到所有的基础解\((x_i, y_i)\)，原方程的正整数解可表示为： \[x + y\sqrt{D} = (x_i + y_i\sqrt{D})(u_1 + v_1\sqrt{D})^n, \quad n \in {0, 1, 2, \dots}\] 如果\(2x^2\)和\(2xy\)中任意一个不能被\(|N|\)整除，则正整数解还包括： \[x + y\sqrt{D} = (x_i - y_i\sqrt{D})(u_1 + v_1\sqrt{D})^n, \quad n \in {1, 2, 3, \dots}\] 可以用递推公式求出基础解对应的解序列： \[\begin{aligned} x’ &= u_1 x + D v_1 y \\ y’ &= v_1 x + u_1 y \end{aligned}\]

用法示例

Project Euler 94: 几乎等边三角形

这题要求的是面积和各边长均为整数的“几乎等边”（两条边相等，第三条边和它们相差不超过1）三角形。

分析：由题意知三边长为\((a,a,b)\)，其中\(a=b\pm1\)。以\(b\)为底边，高\(h = \sqrt{a^2 - \left(\frac{b}{2}\right)^2}\)，面积 \[S = \frac{1}{2} bh = \frac{b\sqrt{4a^2-b^2}}{4} = \frac{b\sqrt{3b^2\pm8b+4}}{4}.\] 要使面积为整数，\(3b^2\pm8b+4\)必须为完全平方数。设 \[3b^2\pm8b+4 = Y^2.\] 两边乘3，得 \[\begin{align*} 9b^2 \pm 24b + 12 &= 3Y^2\\ (3b\pm4)^2 - 4 &= 3Y^2.\\ \end{align*}\] 设\(X=3b\pm4\)，得类Pell方程 \[X^2 - 3Y^2 = 4.\] 注意标准Pell方程\(u^2 - 3v^2 = 1\)的最小正整数解为\((2, 1)\)，因此\(y_\max = \frac{1\cdot\sqrt{4}}{\sqrt{2(2+1)}} < 1\)。这个类Pell方程的基础解只有\((2, 0)\)。 以\((2, 0)\)为起点，用递推式 \[\begin{align*} X’ &= 2X + 3Y\\ Y’ &= X + 2Y \end{align*}\] 求出所有的\((X, Y)\)，由\(X\)反向求出\(b\)。注意：\(b\)必须是正整数，并且还要验证面积S是正整数（经验证，\(Y>0\)时永远满足条件）。

可以用很简短的OCaml代码生成完整的解集：

let solutions =
  Seq.iterate (fun (x, y) -> 2*x + 3*y, x + 2*y) (2, 0)
  |> Seq.filter_map (fun (x, y) ->
      match x mod 3 with
      | 1 when x > 4 -> let b = (x - 4) / 3 in Some (b + 1, b)
      | 2 when y > 0 -> let b = (x + 4) / 3 in Some (b - 1, b)
      | _ -> None)

计算可得如下序列：

a=5 b=6
a=17 b=16
a=65 b=66
a=241 b=240
a=901 b=902
a=3361 b=3360
a=12545 b=12546
...