考虑如下形式的二次丢番图方程：
\[ x^2-Dy^2=1 \]
例如，当\(D=13\)，该方程的最小正整数解是\(649^2-13\times180^2=1\)。不难发现，当\(D\)是完全平方数时，这个方程没有正整数解。对于\(D=\{2,3,5,6,7\}\)，我们可以找到如下最小正整数解：
\[ \begin{aligned}3^2-2\times2^2=1\\2^2-3\times1^2=1\\9^2-5\times4^2=1\\5^2-6\times2^2=1\\8^2-7\times3^2=1\end{aligned} \]
因此对于\(D\le7\)，最小正整数解的中\(x\)的最大值在\(D=5\)时得到。求对于\(D\le1000\)最小正整数解中\(x\)的最大值在\(D\)等于多少时得到？
分析：这道题是到目前为止第一道难度系数为25%的题目，它的难度主要在于题目背后的数学理论，一旦了解了相关理论，这道题的解法还是比较容易理解的。对于熟悉数论的同学，应该比较容易发现题目中所给的方程是一类特殊的丢番图方程，我们一般称之为佩尔方程。关于佩尔方程，数论中已经有非常深入彻底的研究，对其解的存在性以及如何求解都已经有完善的理论。感兴趣的同学可以去翻阅一下相关的数论书籍，我这里推荐希尔弗曼的《数论概论》这本书，这是一本相对入门的数论书籍，其中的第三十章、第三十二章和第四十章专门讨论了佩尔方程，我这里求解佩尔方程的方法也来源于这本书。
佩尔方程的求解方法和六十四与六十五题中涉及的连分数有非常紧密的联系，这里我们也可以看出来欧拉工程安排这种题目顺序的深意。在第六十四题中，我们知道了如何计算非完全平方数平方根的连分数表示以及计算连分数表示的循环周期；在六十五题中，我们学会了如何计算非完全平方数平方根的渐进收敛项，知道可以用连分数的收敛项来对无理数进行近似估计。这些内容是我们在求解佩尔方程中需要用到的，欧拉工程聘用 两道题目为我们预先做了知识储备，这样在做第六十六题的时候不会有过于陡峭的学习曲线。下面我将简述使用平方根的连分数表示及收敛项求解佩尔方程的方法：
对于形如\(x^2-Dy^2=1\)有佩尔方程，当\(D>0\)且为非完全平方数时有无穷多个正整数解，我们这里关心的是最小正整数解，求解步骤如下：
经过上面的步骤就可以求出佩尔方程的最小正整数解，然后只需要在\(D\le1000\)的范围内计算各个方程的解，找到其中最大的\(x\)值，然后求出对应的\(D\)，即为题目所求。事实表明，对于某些\(D\)值，最小正整数解可以变得非常大。如在\(D\le100\)时，最大的最小正整数解在\(D=61\)时取得，它是一个十位数；在\(D\le10000\)时，最大的最小正整数解在\(D=9949\)时取得，它的位数有二百一十二位；在\(D\le100000\)时，最大的最小正整数解在\(D=92821\)时取得，它的位数有七百二十四位！另外，这个网站给出了一个佩尔方程的在线求解器，大家可以代入上面提到的这些值，看看需要多久才能计算出最终结果。
# time cost = 17.2 ms ± 50.2 μs
def pell_equation(D):
a0 = int(D**0.5)
a1 = int(2*a0/(D-a0**2))
seq,P,Q = [a0,a1],[a0,a1*a0+1],[1,a1]
m,d = a0,D-a0**2
while seq[-1] != 2*a0:
m = d*seq[-1]-m
d = (D-m**2)//d
seq.append(int((a0+m)/d))
P.append(seq[-1]*P[-1]+P[-2])
Q.append(seq[-1]*Q[-1]+Q[-2])
if (len(seq)-1) % 2 == 0:
return (P[-2],Q[-2])
else:
return (P[-2]**2+Q[-2]**2*D,2*P[-2]*Q[-2])
def main(N=1000):
dt = {}
for i in range(8,N+1):
if int(i**0.5)**2 == i:
continue
else:
dt[i] = pell_equation(i)[0]
return max(dt,key=dt.get)Project Euler 66: Diophantine equation
原文：https://www.cnblogs.com/metaquant/p/11944011.html