Pell 方程

前置知识：连分数、二次域

引入

本文讨论（广义）Pell 方程的求解。广义 Pell 方程是指关于 \(x\) 和 \(y\) 的不定方程

其中，\(D\) 是正整数且不是完全平方数¹，\(N\) 是非零整数。狭义的 Pell 方程特指 \(N=1\) 或 \(N=\pm 1\) 的特殊情形，有时也包括 \(N=\pm 4\) 的情形。广义 Pell 方程与在实二次整数环内寻找范数为 \(N\) 的二次整数紧密相关，而常常称作（狭义）Pell 方程的这些情形可以看作是寻找实二次整数环内的单位数。

当本文提及 Pell 方程时，特指 \(N=1\) 的情形。相应地，\(N=-1\) 的情形称为负 Pell 方程²（negative Pell's equation）。

解的结构

广义 Pell 方程的整数解 \((x,y)\) 和二次整数 \(x+y\sqrt{D}\) 联系密切，因此文献中 Pell 方程的解常常写作 \(x+y\sqrt{D}\) 的形式。因为二次整数的范数

所以，广义 Pell 方程大致相当于在求解范数为 \(N\) 的二次整数。但是，两者确实有细微的区别。当 \(x\) 和 \(y\) 都是整数时，\(x+y\sqrt{D}\) 一定是二次整数；反过来，二次整数未必要求 \(x\) 和 \(y\) 都是整数——在 \(D\equiv 1\pmod 4\) 的情形，\(x\) 和 \(y\) 还可以同时是半整数³。

这个区别在求解基本单位数时格外重要。因为二次整环中的单位数是指范数为 \(\pm 1\) 的二次整数。对于 \(D\equiv 2,3\pmod 4\)，要找到这样的单位数，只需要求解广义 Pell 方程在 \(N=\pm 1\) 的情形即可；但是对于 \(D\equiv 1\pmod 4\)，还需要考虑 \(N=\pm 4\) 的情形。下文 会讨论单位数的求解方法。

要理解广义 Pell 方程解的结构，需要从 Brahmagupta 恒等式 入手：

它相当于二次整数的范数保持乘法，即

利用这一恒等式，可以利用方程 \(x^2-Dy^2=N_1\) 和方程 \(x^2-Dy^2=N_2\) 的整数解复合出方程 \(x^2-Dy^2=N_1N_2\) 的整数解。当然，从二次整数的角度看，解的复合就是二次整数的乘法，这就体现了将 Pell 方程的解记作二次整数形式的方便之处。特别地，取 \(N_1=N\) 和 \(N_2=1\) 就可以发现，如果已知 \(x^2-Dy^2=N\) 的一组解和相应的 Pell 方程 \(x^2-Dy^2=1\) 的全体解，就可以得到 \(x^2-Dy^2=N\) 更多的解。当然，这种方法未必能够生成全部的解。但是，这至少说明理解 Pell 方程的解的结构对理解广义 Pell 方程的解的结构有重要作用。

Pell 方程

方程 \(x^2-Dy^2=1\) 的几何意义是实轴在 \(x\) 轴、虚轴在 \(y\) 轴的双曲线。双曲线上的每个点都唯一对应了 \(x+y\sqrt{D}\) 的一个非零取值：双曲线的左支对应着 \(x+y\sqrt{D}\) 的负值，右支则对应正值。而且，在每一支上，双曲线自下而上对应的 \(x+y\sqrt{D}\) 的取值是严格递增的。二次整数的取值给 Pell 方程的解赋予了自然的顺序。

双曲线同时关于 \(x\) 轴和 \(y\) 轴对称，因此讨论 Pell 方程的解只需要考虑在第一象限内的那一段即可，其余解可以通过对称性获得。这相当于只考虑 \(x+y\sqrt{D}>1\) 的解。如果方程除了 \((\pm 1,0)\) 之外还存在不平凡的解，那么一定在第一象限内存在 \(x+y\sqrt{D}\) 取值最小的解 \((x_1,y_1)\)，这也是第一象限内（不含坐标轴）横纵坐标都最小的整点，它称为 Pell 方程的基本解（fundamental solution）⁴。根据前文的讨论，满足 \(x_k+y_k\sqrt{D}=(x_1+y_1\sqrt{D})^k\) 的整数对 \((x_k,y_k)\) 都是 Pell 方程的解，而且都在第一象限。反过来，这也的确是 Pell 方程在第一象限内的全部解。再利用对称性，就可以得到如下结论：

前文的讨论只是假设了基本解的存在。现在要说明的是，Pell 方程总是存在非平凡的解。

当然，本节提供的是非构造性的证明，在下文讨论 Pell 方程的解法时，会直接利用连分数的渐近分数构造出 Pell 方程的解，因而提供了 Pell 方程存在非平凡解的另一种证明。另外，尽管本节得到的 Pell 方程解的结构与实二次整数环的单位数的结构是一致的，但是对于 \(D\equiv 1\pmod 4\) 的情形，本节并没有完全解决相应的二次整数环的单位数的结构问题，下文将进一步讨论。

广义 Pell 方程

广义 Pell 方程 \(x^2-Dy^2=N\) 的图像同样是平面上的双曲线，同样以 \(x\) 轴和 \(y\) 轴为对称轴。前文已经指出，方程 \(x^2-Dy^2=N\) 的部分解可能只相差一个 Pell 方程解的因子，这意味着可以将方程 \(x^2-Dy^2=N\) 的解划分为等价类。对于方程 \(x^2-Dy^2=N\) 的两个解 \((x_1,y_1)\) 和 \((x_2,y_2)\)，如果存在 Pell 方程的解 \((u,v)\) 使得 \(x_2+y_2\sqrt{D}=(x_1+y_1\sqrt{D})(u+v\sqrt{D})\) 成立，那么称解 \((x_1,y_1)\) 和 \((x_2,y_2)\) 等价。两个解等价的充分必要条件是

因为 Pell 方程的解相对容易求出，一个自然的想法是在上述的每个等价类中各求出一个解。只要知道这些解，就可以利用相应的 Pell 方程的解得到所要求的广义 Pell 方程的全部解。在广义 Pell 方程的解的等价类中，由于对称性，每个等价类都存在纵坐标 \(y\) 非负但是尽可能小的解：如果这样的解唯一，它就称为该等价类的基本解；否则，该等价类必然有两个 \(y\) 非负且最小的解，而且它们关于 \(y\) 轴对称，此时选择 \(x>0\) 的那个作为基本解。由此，求解广义 Pell 方程 \(x^2-Dy^2=N\)，就相当于求解它的基本解集 \(U\)。设它对应的 Pell 方程的基本解是 \((r,s)\)，则广义 Pell 方程的全部解的集合是

广义 Pell 方程的基本解必然是有限的。这是因为从上面的通解表达式可知，绝对值 \(|u+v\sqrt{D}|\) 必然位于 \(r-s\sqrt{D}\) 和 \(r+s\sqrt{D}\) 之间。文末的参考文献中提供了关于基本解的坐标的范围的更严格的估计。当然，与 Pell 方程的情形不同，广义 Pell 方程可能没有解。

利用广义 Pell 方程的一个解 \((u,v)\) 和 Pell 方程的基本解 \((r,s)\) 得到同一个等价类中的全部解的方法，除了利用解的复合之外，还可以利用如下递推关系

其中，\(x_k+y_k\sqrt{D}=(r+s\sqrt{D})^k(u+v\sqrt{D})\)。这是因为 \(x_n\) 和 \(y_n\) 都可以对某一对实数 \((A,B)\) 写成 \(A(r+s\sqrt{D})^k+B(r-s\sqrt{D})^k\) 的形式，而根据 Vieta 定理，\(r\pm s\sqrt{D}\) 是方程 \(x^2-2rx+1=0\) 的两个实根，进而 \(x_n\) 和 \(y_n\) 都满足上述的二阶常系数递推关系。相较于解的复合，该递推公式有更少的乘法次数。

求解方法

Pell 方程和广义 Pell 方程的求解都可以基于连分数进行。

PQa 算法

本文讨论的算法都基于 PQa 算法，它可以用于求出特定的二次无理数的连分数展开。

设整数 \(P_0,Q_0,D\) 满足 \(Q_0\neq 0\)，\(D>0\) 且不是完全平方数，以及 \(P_0^2\equiv D\pmod{Q_0}\)。那么，二次无理数

的连分数展开 \([a_0,a_1,\cdots]\) 可以通过如下 递推公式 求得：

进而，\(\omega\) 的第 \(k\) 个渐近分数的分子和分母 \(A_k\) 和 \(B_k\) 由如下 递推公式 给出：

且 \(A_{-1} = 1\)，\(A_{-2}=0\)，\(B_{-1}=0\)，\(B_{-2}=1\)。

这些公式的正确性已经在连分数一文中得到证明。那里也说明了，因为二次无理数是 循环连分数，所以，三元组 \((P_k,Q_k,a_k)\) 最终将进入循环，算法总可以在有限步内终止。不妨设循环节的最小长度是 \(\ell\)，且循环的最早的起始位置是 \(k_0\)，则二次无理数的连分数展开可以写作

利用 PQa 算法解决 Pell 方程，需要建立如下结论：

这个结论提供了一种寻找方程 \(x^2-Dy^2=N\) 的解的方法。如果合理地选择 \(Q_0>0\) 并选择 \(P_0\) 为同余方程 \(P_0^2\equiv D\pmod{Q_0}\) 的解，那么通过对 \((P_0+\sqrt{D})/Q_0\) 执行 PQa 算法，直到找到 \((-1)^kQ_0Q_{k}=N\)，此时 \((G_{k-1},B_{k-1})\) 就成为原方程的一组解。而且，如果 \(Q_k=\pm 1\)，那么这样得到的解一定是本原解，也就是说 \(G_{k-1}\) 和 \(B_{k-1}\) 一定是互素的。

这个思想是解决 Pell 方程和广义 Pell 方程的核心。理解了这一思想后，下面就着手处理算法的一些细节，并证明所有的解都可以通过该方式得到。

Pell 方程

要解决 Pell 方程 \(x^2-Dy^2=1\)，只需要对 \((P_0,Q_0,D)=(0,1,D)\) 运行 PQa 算法，直到出现 \((-1)^kQ_k=1\)，此时 \((A_{k-1},B_{k-1})\) 就是 Pell 方程的一组解（因为 \(G_{k-1}\) 此时就是 \(A_{k-1}\))。当然，对于 Pell 方程，对该过程可以进行更精确的描述。

首先，解一定出现在循环节的末尾处。上述过程相当于对 \(\sqrt{D}\) 做连分数展开。对此，已经有 结论：

此处，循环节长度为 \(\ell\)，且起始位置是第 \(1\) 项（下标从 \(0\) 开始）。而且，它的第 \(\ell\) 项余项等于 \(\lfloor\sqrt{D}\rfloor+\sqrt{D}\)，这说明 \(Q_{\ell}=1\)。因此，如果 \(\ell\) 是偶数，那么 \((A_{\ell-1},B_{\ell-1})\) 就是 Pell 方程的一组非平凡解；如果 \(\ell\) 是奇数，那么 \((A_{2\ell-1},B_{2\ell-1})\) 就是 Pell 方程的一组非平凡解。

接下来要说明，刚刚得到的这组解一定是基础解。这个结论基于两点理由：第一，Pell 方程的所有正整数解 \((x,y)\) 对应的分数 \(x/y\) 都出现在 \(\sqrt{D}\) 的渐近分数中，这就保证了 \((x,y)\) 必定是 PQa 算法过程中的某个 \((A_k,B_k)\)；第二，除了循环节末尾，不会再出现其他位置有 \(Q_k=1\)，因为 \(A_k\) 和 \(B_k\) 的递推关系保证了它们的大小随着下标增加而增加，所以最小的正整数解（即基础解）必然出现在刚刚指明的位置。这两点理由可以分别从如下的两条定理得出：

综合本节的讨论可知，只要对 \(\sqrt{D}\) 做连分数展开，亦即以 \((P_0,Q_0,D)=(0,1,D)\) 为起点做 PQa 算法，当首次得到 \(Q_\ell=1\) 时，就到达了第一个循环节的末尾。此时，如果 \(\ell\) 是偶数，那么 \((A_{\ell-1},B_{\ell-1})\) 就是 Pell 方程的基本解；否则，\((A_{2\ell-1},B_{2\ell-1})\) 是 Pell 方程的基本解。对于循环节长度 \(\ell\) 为奇数的情形，并不需要继续 PQa 算法到两倍的循环节处，马上就会说明 \(A_{2\ell-1}+B_{2\ell-1}\sqrt{D}=(A_{\ell-1}+B_{\ell-1}\sqrt{D})^2\)，因而可以直接从 \((A_{\ell-1},B_{\ell-1})\) 直接计算出 Pell 方程的基本解。Pell 方程所有其他解都可以通过 Pell 方程的基本解计算。

负 Pell 方程

根据上一节的讨论可知，负 Pell 方程的解也必然对应于 \(\sqrt{D}\) 的渐近分数，而且只能出现在 \((-1)^kQ_k=-1\) 处。这只能出现在循环节的末尾。因此，负 Pell 方程有解，当且仅当循环节长度 \(\ell\) 是奇数。当解存在时，\((A_{\ell-1},B_{\ell-1})\) 就是负 Pell 方程的基本解。它的求解方法和上一节是一致的。

利用前文对于 Pell 方程解的结构的证明相仿的思路，可以证明如下结论：

因为 \((A_{\ell-1},B_{\ell-1})\) 是负 Pell 方程的最小正整数解，而 \(x^2-Dy^2=\pm 1\) 的所有正整数解都出现在集合

中，又因为这些正整数解必然对应 \(\sqrt{D}\) 的在循环节末尾（前一位）处的渐近分数，而且渐近分数的分子和分母是严格单调递增的，所以对所有 \(k\in\mathbf N_+\) 总是有

在所有这些正整数解中，\(k\) 为奇数时就是负 Pell 方程的解，\(k\) 为偶数时就是 Pell 方程的解，两者交替出现。

判断负 Pell 方程是否有解，需要计算 \(\sqrt{D}\) 连分数展开的循环节的长度，这并不容易计算，因此希望能够找到更简单的判断方法。但是，目前并没有条件简明、容易计算的判断方法⁵。此处仅仅提供一个简单的结论。

如果 \(D\) 是合数，那么不含有 \(4k+3\) 型素因子且没有平方因子也不能保证方程 \(x^2-Dy^2=-1\) 有解，例如 \(x^2-34y^2=-1\) 就没有解。

范数为 ±4 的情形

接下来讨论方程 \(x^2-Dy^2=\pm 4\) 的解。此时，解的性态取决于 \(D\bmod 4\) 的大小。

有些情形是容易的。如果 \(D\equiv 0\pmod 4\)，那么 \(x\) 是偶数，因而 \((x/2,y)\) 是方程 \(u^2-(D/4)v^2=\pm 1\) 的解。其余的情形，必然有 \(x,y\) 同时是奇数或者同时是偶数。如果 \(x,y\) 同时是奇数，方程两侧对 \(4\) 取模就得到 \(D\equiv 1\pmod 4\)。所以，如果 \(D\equiv 2,3\pmod 4\)，那么 \(x,y\) 只能同时是偶数，因而 \((x/2,y/2)\) 是方程 \(u^2-Dv^2=\pm 1\) 的解。因此，除了 \(D\equiv 1\pmod 4\) 的情形，方程 \(x^2-Dy^2=\pm 4\) 的解都可以通过相应的（负）Pell 方程的解得到。

现在讨论 \(D\equiv 1\pmod 4\) 的情形，它不能简单地转化为已经解决的情形。为了找到基本解，可以对 \((P_0,Q_0,D)=(1,2,D)\) 应用 PQa 算法，当首次得到 \(Q_\ell=2\) 时，就到达了第一个循环节的末尾。如果循环节长度 \(l\) 是偶数，那么 \((G_{\ell-1},B_{\ell-1})\) 就是方程 \(x^2-Dy^2=4\) 的基本解；否则，\((G_{\ell-1},B_{\ell-1})\) 就是方程 \(x^2-Dy^2=-4\) 的基本解。从 \((G_{\ell-1},B_{\ell-1})\) 出发，可以得到方程 \(x^2-Dy^2=\pm 4\) 的所有解：

如果循环节长度 \(\ell\) 是偶数，所有这些都是方程 \(x^2-Dy^2=4\) 的解；否则，当 \(k\) 是奇数时，\((x,y)\) 是方程 \(x^2-Dy^2=-4\) 的解，而当 \(k\) 是偶数时，\((x,y)\) 是方程 \(x^2-Dy^2=4\) 的解。

这个算法的正确性依赖于如下事实：

综合这些事实，重复前文几节的论述，就可以说明用于解决方程 \(x^2-Dy^2=\pm 4\) 的上述算法的正确性。这些结果说明，方程 \(x^2-Dy^2=\pm 4\) 具有和方程 \(x^2-Dy^2=\pm 1\) 类似的简单的解的结构：它的所有解都可以通过其最小正整数解表示出来，而无需求解其它方程。

其实，方程 \(x^2-Dy^2=\pm 1\) 的所有解都可以在方程 \(x^2-Dy^2=\pm 4\) 的解中找到，因而从这个角度看，方程 \(x^2-Dy^2=\pm 4\) 更为基础。显然，\((x,y)\) 是方程 \(x^2-Dy^2=\pm 1\) 的解，当且仅当 \((2x,2y)\) 是方程 \(x^2-Dy^2=\pm 4\) 的解。前文的分析指出，当 \(D\equiv 2,3\pmod 4\) 时，方程 \(x^2-Dy^2=\pm 4\) 的所有解都一定同为偶数，因而对应于 \(x^2-Dy^2=\pm 1\) 的解。

当 \(D\equiv 0\pmod 4\) 时，方程 \(x^2-Dy^2=\pm 4\) 的解 \((x,y)\) 中，\(x\) 一定是偶数，但是 \(y\) 可能是奇数。如果在方程 \(x^2-Dy^2=\pm 4\) 的最小正整数解 \((x_1,y_1)\) 中，\(y_1\) 是偶数，那么在所有解中 \(y\) 也一定是偶数，此时这些整数解和方程 \(x^2-Dy^2=\pm 1\) 的整数解一一对应；但是，如果在最小整数解 \((x_1,y_1)\) 中，\(y_1\) 是奇数，那么 \(y_k\) 的奇偶性将和 \(k\) 一致，交替变化，因而只有当 \(k\) 是偶数时，才对应于方程 \(x^2-Dy^2=\pm 1\) 的解。如果 \(x^2-Dy^2=\pm 4\) 的最小正整数解中 \(y_1\) 是奇数而且 \(x_1+y_1\sqrt{D}\) 的范数是 \(-4\)，那么，对于这样的 \(D\)，\(x^2-Dy^2=-4\) 有解，但是 \(x^2-Dy^2=-1\) 无解。

当 \(D\equiv 1\pmod 4\) 时，方程 \(x^2-Dy^2=\pm 4\) 的解 \((x,y)\) 可能同时是奇数，也可能同时是偶数。如果最小正整数解 \((x_1,y_1)\) 已经同时是偶数，那么它的所有整数解也一定同时是偶数，所以总是对应于方程 \(x^2-Dy^2=\pm 1\) 的整数解。如果最小正整数解 \((x_1,y_1)\) 同时是奇数，那么有如下结论：

也就是说，方程 \(x^2-Dy^2=\pm 4\) 的每三个解中就有一个同时是偶数，它对应着 \(x^2-Dy^2=\pm 1\) 的整数解。这也说明，对于 \(D\equiv 1\pmod 4\)，方程 \(x^2-Dy^2=-4\) 有解，当且仅当方程 \(x^2-Dy^2=-1\) 有解。

到目前为止的讨论，已经足够计算实二次整数环的基本单位数。设 \(D\) 是正整数且不含平方因子。对于 \(D\equiv 2,3\pmod 4\) 的情形，只需要求出 \(x^2-Dy^2=\pm 1\) 的最小正整数解；而对于 \(D\equiv 1\pmod 4\) 的情形，只需要求出 \(x^2-Dy^2=\pm 4\) 的最小正整数解。当得到最小正整数解 \((x,y)\) 时，对于 \(D\equiv 2,3\pmod 4\)，基本单位数就是 \(\pm x\pm y\sqrt{D}\)；对于 \(D\equiv 1\pmod 4\)，基本单位数就是 \(\dfrac{\pm x\pm y\sqrt{D}}{2}\)。

一般情形

最后，讨论广义 Pell 方程的解法。

对于 \(|N|<\sqrt{D}\) 的情形有一个简单的解法。前文的结论说明方程 \(x^2-Dy^2=N\) 的解 \((x,y)\) 一定满足 \(\dfrac{x}{y}\) 等于 \(\sqrt{D}\) 的某个渐近分数。而且，根据前文讨论的解的结构可知，每个基础解 \((x,y)\) 都满足 \(x+y\sqrt{D}\) 小于等于相应的 Pell 方程 \(x^2-Dy^2=1\) 的基础解 \(x_1+y_1\sqrt{D}\)。利用 PQa 算法中分母序列 \(B_k\) 的单调性可知，广义 Pell 方程的这些基础解一定出现在相应的 Pell 方程的基础解出现之前。由此，只需要对 \((P_0,Q_0,D)=(0,1,D)\) 运行 PQa 算法，直到 \(Q_{\ell'}=1\) 且 \(\ell'\) 为偶数时为止，过程中对出现的每个 \((A_k,B_k)\) 检验是否存在整数 \(f\) 使得

成立，如果成立，则记录 \((fA_{k},fB_{k})\) 是一个最小正整数解。这个过程记录的所有 \((fA_k,fB_k)\) 就是方程 \(x^2-Dy^2=N\) 的全部最小正整数解。利用 \((A_{\ell'-1},B_{\ell'-1})\)，即相应的 Pell 方程的基本解，就可以根据求出的这些最小正整数解，生成广义 Pell 方程的所有解。注意，取决于循环节长度 \(\ell\) 是偶数还是奇数，上述的 \(\ell'\) 可能是 \(\ell\) 或是 \(2\ell\)。

对于更为一般的 \(N\) 的情形，上述方法不再适用。首先，枚举 \(N\) 的所有平方因子 \(f^2\)，设 \(m=N/f^2\)，并枚举同余方程 \(z^2\equiv D\pmod{|m|}\) 的所有满足 \(-|m|/2<z \le |m|/2\) 的解 \(z\)。然后，对 \((P_0,Q_0,D)=(z,|m|,D)\) 运行 PQa 算法，直到 \(Q_k=\pm 1\) 或已经结束了一个循环节。在第二种情形，那么与该组 \((f,z)\) 相关的方程的解并不存在。在第一种情形，需要进一步判断 \((-1)^kQ_k=N/|N|\) 与否。如果符号一致，那么 \((fG_{k-1},fB_{k-1})\) 就是方程 \(x^2-Dy^2=N\) 的解。否则，它是方程 \(x^2-Dy^2=-N\) 的解，而且当且仅当相应的负 Pell 方程的解存在时，才可以通过复合它与相应的负 Pell 方程的基本解来得到方程 \(x^2-Dy^2=N\) 的解。当完成对所有组 \((f,z)\) 的遍历之后，就可以得到方程 \(x^2-Dy^2=N\) 在每个解的等价类中各恰好一个解，且该解为该等价类中的基本解或最小正整数解。利用它们和相应的 Pell 方程的基本解，可以生成该方程的所有整数解。这一算法称为 Lagrange–Matthews–Mollin 算法。

该算法的正确性由如下定理保证：

该定理保证了方程的所有正解都存在于相应的二次无理数的渐近分数中。因为利用 PQa 算法计算渐近分数时，只要进入循环节，就一定能保证渐近分数总是正的。所以，只要枚举所有定理条件所允许的二次无理数，计算它的渐近分数直到一个循环节内，就能够找到一个解。因为在同一个二次无理数的渐近分数中出现的两个解，一定是等价的，所以只要得到第一个满足 \((-1)^kQ_k=N/|N|\) 的解，就可以停止继续后面的计算。与前面的所有算法都不同的是，此处满足条件的 \(k\) 可能出现在尚未进入循环节时。

习题

• LOJ 6687.「Project Euler 66」解方程
• SPOJ EQU2 - Yet Another Equation
• SPOJ PELL2 - Pell (Mid pelling)
• UVa 12909. Numeric Center
• UVa 10241. Semi-triangular and also Square

参考文献与注释

• Pell's equation - Wikipedia
• John P. Robertson - Solving the generalized Pell equation \(x^2-Dy^2=N\)
• Keith Matthews - The Diophantine Equation \(x^2-Dy^2=N\),\(D>0\)
• Existence of Solution to Pell’s Equation - Suryateja Gavva's Blog
• Calculating the simple continued fraction of a quadratic irrational - Number Theory Web（PQa 算法）
• Solving the diophantine equation x2–Dy2 = N, D > 0 and not a perfect square, N ≠ 0 - Number Theory Web（Lagrange–Matthews–Mollin 算法）

本页面的全部内容在 小熊老师 - 莆田青少年编程俱乐部 0594codes.cn 协议之条款下提供，附加条款亦可能应用