Pell方程

考虑如下形式的二次丢番图方程（即Pell方程）的正整数解，其中D是自然数： $$x^2-Dy^2=1.$$ 显然，D是完全平方数时不存在正整数解。

D不是完全平方数正整数解总是存在，例如$D=13$时有$649^2 – 13\times 180^2 = 1$，这也$D=13$时是x值最小的一组解。

$D=2,3,5,6,7$时，x值最小的解分别是 \begin{align*} 3^2 - 2\times 2^2 &= 1 \\ 2^2 - 3\times 1^2 &= 1 \\ {\color{red}9}^2 - 5\times 4^2 &= 1 \\ 5^2 - 6\times 2^2 &= 1 \\ 8^2 - 7\times 3^2 &= 1 \\ \end{align*} 因此对不超过7的D而言，$D=5$时x的最小值最大。

考虑不超过1000的D，问D取何值时x的最小值最大。

本题难度:

解答

参考Pell方程的解法，x值最小的正整数即基本解，可由$\sqrt D$的连分数展开得到，结果是$D=661$。

注：此时$x=16421658242965910275055840472270471049$，$y=638728478116949861246791167518480580$。

import math

sq=[i*i for i in range(33)]
x=y=d=0
for n in range(2,1001):
    if n not in sq:
        s=int(math.sqrt(n))
        r,k,a=s,1,[s]
        while k!=1 or len(a)==1:
            k=(n-r*r)/k
            r=(s+r)/k*k-r
            a.append((s+r)/k)
        if len(a)%2:
            a=a[:-1]
        else:
            a=a+a[1:-1]
        p,q,i=[a[0],a[0]*a[1]+1],[1,a[1]],2
        while i < len(a):
            p.append(a[i]*p[-1]+p[-2])
            q.append(a[i]*q[-1]+q[-2])
            i=i+1    
        if p[-1]>x:
            x,d,y=p[-1],n,q[-1]
print x,d,y