Project Euler Problem 066

Statement

onsider quadratic Diophantine equations of the form:

x2 – Dy2 = 1

For example, when D=13, the minimal solution in x is 6492 – 131802 = 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:

32 – 2*22 = 1
22 – 3*12 = 1
92 – 5*42 = 1
52 – 6*22 = 1
82 – 7*32 = 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.

Solution

The solution came out after reading and understanding this documents:
Wikipedia Pell's Equation
Wikipedia Convergent in continued fractions
Wikipedia Fundamental Recurrence Formulas

It's a basic implementation of all stated in those three articles.

def continued_fraction(n):
    m = 0
    d = 1
    a = a_0 = int(n ** 0.5)
    triple_set = list()
    while (m, d, a) not in triple_set:
        triple_set.append((m, d, a))
        tmp_m = d * a - m
        tmp_d = (n - tmp_m ** 2) // d
        tmp_a = int((a_0 + tmp_m) / tmp_d)
        a = tmp_a
        d = tmp_d
        m = tmp_m
    triple_set.append((m, d, a))
    return triple_set
 
def min_x(n):
    resolution = continued_fraction(n)
    restart_index = resolution.index(resolution[-1])
    index = min(2, len(resolution) - 1)
    stored_x = []
    stored_y = []
    x = resolution[0][2]
    y = 1
    if x ** 2 - n * (y ** 2) == 1:
        return x
    stored_x.append(x)
    stored_y.append(y)
    x = resolution[0][2] * resolution[1][2] + 1
    y = resolution[1][2]
    while x ** 2 - n * (y ** 2) != 1:
        stored_x.append(x)
        stored_y.append(y)
        x = resolution[index][2] * x + 1 * stored_x[-2]
        y = resolution[index][2] * y + 1 * stored_y[-2]
        index += 1
        if index == len(resolution):
            index = restart_index + 1
    return x
 
if __name__ == '__main__':
    max_x = 0
    result = 0
    for d in filter(lambda d: int(d ** 0.5) ** 2 != d, range(1, 1001)):
        tmp = min_x(d)
        if tmp > max_x:
            max_x = tmp
            result = d
    print("The result is:", result)

The Python file is available for download here.

Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License