Problem 66

Problem 66

26 March 2004

Consider quadratic Diophantine equations of the form:

x² – Dy² = 1

For example, when D=13, the minimal solution in x is 649² – 13180² = 1.

It can be assumed that there are no solutions in positive integers when D is square.

By finding minimal solutions in x for D = {2, 3, 5, 6, 7}, we obtain the following:

3² – 22² = 1
2² – 31² = 1
9² – 54² = 1
5² – 62² = 1
8² – 73² = 1

Hence, by considering minimal solutions in x for D 7, the largest x is obtained when D=5.

Find the value of D 1000 in minimal solutions of x for which the largest value of x is obtained.

x² – Dy² = 1这个方程叫做pell方程，对于非完全平方数D, X, Y有无数个解

求解方法 一般用的是连分数法

求出SQRT(D)的连分数[a₀:(a₁, a₂, a₃ ... a_n)...],其中a_n=2*a₁

然后再计算

根据连分数循环节位数n的奇偶性分两种情况

若n为偶数,则 x= p_n-1 , y= q_n-1 ; n为奇数时x= p_2n-1 , y= q_2n-1

p₀= a₀, p₁= a₀* a₁+1 , p_n= a_n*p_n-1+p_n-2
q₀=1,   q₁=a₁,        q_n=a_n*q_n-1+q_n-2

注:以上说的解均为最小正解

def getSqrtFraction(num):
    result = []
    z     = int(num**0.5)
    result = [z]
    if z*z == num:
        return result

    d     = {}#用于保存中间已计算变量
    fm    = (num, z)
    fz    = 1
    while 1:
        k = [fm, fz]
        if d.has_key(str(k)):
            return result
        else:
            d[str(k)] = True
            t = (num - z**2)/fz
            result.append(int((num**0.5+z)/t))
            z = (int((num**0.5+z)/t))*t - z
            fm = (num, z)
            fz = t


def fun66():
    result, rnum = 0, 0
    for i in range(2, 1001):
        fc = getSqrtFraction(i)

        k  = len(fc)
        if k  == 0:
            continue
        if k == 1:
            continue
        p0 = fc[0]
        p1 = fc[0]*fc[1]+1
    
        n  = k - 1
        if n%2 == 0:#循环链为偶数
            t  = 2
            pn = p1
            while t < n:
                pn = fc[t]*p1 + p0
                p0, p1 = p1, pn
                t += 1
            if result < pn:
                result = pn
        else:#循环链长度为奇数
            t  = 2
            pn = p1
            while t < 2*n:
                if t <= n:
                    pn = fc[t]*p1 + p0
                else:
                    pn = fc[t%n]*p1 + p0
                p0, p1 = p1, pn
                t += 1
            if result < pn:
                rnum   = i
                result = pn
    return result, rnum

getSqrtFraction(num)用于获取num的连分数
answer is 661, x=16421658242965910275055840472270471049
time:0.125