Project Euler Problem 66 Solution with python

Diophantine equation¶

Consider quadratic Diophantine equations of the form:

$$ x^{2} – Dy^{2} = 1 $$

For example, when $D = 13$, the minimal solution in $ x $ is $ 649^{2} – 13 \times 180^{2} = 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^{2} – 2×2^{2} = 1 $$ $$ 2^{2} – 3×1^{2} = 1 $$ $$ 9^{2} – 5×4^{2} = 1 $$ $$ 5^{2} – 6×2^{2} = 1 $$ $$ 8^{2} – 7×3^{2} = 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.

If you have already solved Problem 64 on Project euler, then assume that this one is an extension to Problem 64.

If you have not solved Problem 64, try solving that problem first. Of course I have written a post on that.

The solution is based on two approaches: Chakravala method(Page No:3) and Continued Fraction method(Page No: 18). And to find continued fractions I have used this source: Wikipedia - Methods of Computing square roots.

As there is a lot of stuff in the given sources, the screen grabs of the required stuff is as follows:

Theorem statement to solve pells equation

Bhaskara's Lemma to solve negative pell's equation

Algorithm to find the continued fractions for square root of a number

There are four things that we will have to understand from the sources:

1. If the period of the continued fraction is even then we are finding the solution for the positive pells equation - $ x^{2} - Dy^{2} = 1 $
2. If the period of the continued fraction is odd then we are finding the solution for the negative pells equation - $ x^{2} - Dy^{2} = -1 $
3. What ever might be the case, we don't need the last value of the convergent list. i.e. in our case $ \sqrt{19} = [4;2,1,3,1,2] $ instead of $ [4;2,1,3,1,2,8] $
4. If we are solving the negative pells equation, then we will use Baskara's Lemma to get the required answer.

Now that you have understood the important points, have a look at the program and you are good to go.

Program¶

In [ ]:

# http://radiusofcircle.blogspot.com

#import time
import time

# square root function
from math import sqrt

# time at the start of program execution
start = time.time()

# function to calculate the continued fraction
def cf(n):
    mn = 0.0
    dn = 1.0
    a0 = int(sqrt(n))
    an = int(sqrt(n))
    convergents = [a0]
    period = 0
    if a0 != sqrt(n):
        while an != 2*a0:
            mn = dn*an - mn
            dn = (n - mn**2)/dn
            an = int((a0 + mn)/dn)
            convergents.append(an)
    return convergents[:-1]

def cf_inv(cf):
    """
    function to calculate the
    simple fraction from the continued
    fraction.
    """
    numerator = 1
    denominator = cf.pop()
    while cf:
        denominator, numerator = denominator*cf.pop() + numerator, denominator
    return denominator, numerator

# variable to store the largest value 
# and the place it occurs
largest = 0, 0

# for loop less than 1000
for i in xrange(1, 1001):
    if i%sqrt(i) != 0:
        continued_fraction = cf(i)
        if len(continued_fraction) % 2 != 0:
            u, v = cf_inv(continued_fraction)
            u, v = 2*u**2+1, 2*u*v
        else:
            u, v = cf_inv(continued_fraction)
        if u > largest[1]:
            largest = i, u

# print the largest value
print largest[0]

# time at the end of program execution
end = time.time()

# total time taken to run the program
print end - start

I have tried to name the variables as meaningful as possible. But a brief explanation is as follows:

The function cf will take a number($ n $) as input and calculate the continued fraction for square root of the number($ \sqrt{n} $).

So the expected output will be [a0, a1, a2, a3, .........., an]. Also from the reference, we know that an = 2a0. We also know that, to solve pell's equation we don't need the last value(a_n) and so cf will return a list with values excluding the last number.

The function cf_inv will take a list(List of continued fraction values) as input and will return the numerator and the denominator if the list is written as a simple fraction.

For example, consider the list: [2, 3, 4, 5]. This list when written as a continued fraction will be as follows: \begin{equation} x = 2 + \cfrac{1}{3 + \cfrac{1}{4 + \cfrac{1}{5} } } \\ \implies x = 2 + \cfrac{1}{3 + \cfrac{5}{21} } \\ \implies x = 2 + \cfrac{21}{68} \\ \implies x = \cfrac{157}{68} \end{equation} Now that you have understood how the function works in the background, we will use tweak the function so that it will print the numerator and denominator for each and every iteration.

In [1]:

def cf_inv(cf):
    numerator = 1
    denominator = cf.pop()
    while cf:
        denominator, numerator = denominator*cf.pop() + numerator, denominator
        print numerator, denominator
    return denominator, numerator

print cf_inv([2, 3, 4, 5])

5 21
21 68
68 157
(157, 68)

It should be noted that, when using the while loop, during the last iteration we will swap the numerator and denominator(Think why I am doing that and comment in the comment box below).

Rest of the code is simple and easy to understand, just a for loop to find the answer.

Output¶

Answer: 661
Execution time: 0.00883007049561 seconds
Output:

Project Euler Problem 66 Terminal output

Project Euler Problem 66 with answer

Summary¶

When I started solving this problem, I started reading papers and articles on chakravala method and tried implementing the same. But for a few numbers I always got a wrong answer. I didn't knew why! After a few trial and errors I again started searching on the web, but this time "Solving Pell's equation" and finally found perfect source to solve the question. Further discussing about the problem, it was not really easy and I personally liked the problem. Still there is a room for improvement. But for now I am satisfied with the speed of execution.

If you have any doubts(not just related to this problem but any general ones are also welcome) or didn't understand anything comment in the comment box below and I will be glad to help you. Please comment in the comment box below if you have found any typo or have a better program.

You can also contact me.

Thank you. Have a nice day😃.