Problem 66 - Investigate the Diophantine equation x2 − Dy2 = 1.

링크

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.

It was hardest problem I had solved from Project Euler ever.

First approach - rather naive one

Algorithm: Check for all integer x >= 2 until the equation hold.

This approach didn't work efficiently because the variable x was only increased by 1.

In case d=61, the minimum possible variable for x is 1766319049. and I couldn't find the number using this algorithm.

Moreover, in case d=661, the smallest x is 16421658242965910275055840472270471049 !! that is the number this problem wants to find and is amazingly enormous number. 

This problem cannot be solved by this kind of algorithms(naive algorithms).

To solve this problem, we must find pell's equation.

Second approach - Pell's equation and Continued Fraction

Let  denote the sequence of convergents to the continued fraction for . Then the pair (x₁,y₁) solving Pell's equation and minimizing x satisfies x₁ =h_i and y₁ = k_i for some i. This pair is called the fundamental solution. Thus, the fundamental solution may be found by performing the continued fraction expansion and testing each successive convergent until a solution to Pell's equation is found.

-from wikipedia

Once we know the theory, it is easy to implement

python

import math maxX, maxD = 0, 0 for d in range(1, 1001): k = int(d**0.5) if k*k == d: continue a, b = k, 1 hn, kn = k, 1 hp, kp = 1, 0 while hn*hn - d*kn*kn != 1: p = (d - a*a) / b q = int(b / (d**0.5 - a)) r = p*q - a hn, hp = q*hn + hp, hn kn, kp = q*kn + kp, kn a, b = r, p if maxX < hn: maxX, maxD = hn, d print maxX, maxD