(24 votes, average: 5.00 out of 5)

## Project Euler 66: Investigate the Diophantine equation x2 − Dy2 = 1.

#### Problem Description

Consider quadratic Diophantine equations of the form:

x2 – Dy2 = 1

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

#### Analysis

Using the chakravala method for solving minimal solutions to Pell’s Equation. A detailed explanation is posted in the comments section.

#### Project Euler 66 Solution

Runs < 0.005 seconds in Python 2.7.

Use this link to get the Project Euler 66 Solution Python 2.7 source.

#### Afterthoughts

• D=92821 for D ≤ 100000 in 15 seconds.

Solution of “Pell’s Equation” by Chakravala Method by Gopal Menon

It was discovered by Brahmagupta that a triple (a, b, k), representing the equation a2 – Nb2 = k, can be composed with the trivial triple (m, 1, m2 – N) to get a new triple (am + Nb, a + bm, k(m2 – N)). Scaling down by k, using Bhaskara’s lemma, we get,

a2 – Nb2 = k ==> {(am + Nb)/k}2 – N{(a + bm)/k}2 = (m2 – N)/k

The following is the algorithm for solving “Pell’s equation”, x2 – Ny2 = 1, as given by the Indian mathematicians who were the first to solve this equation:

1. For the given value of N (which is not a perfect square) replace x by a, y by b and 1 by k and solve the equation a2 – Nb2 = k. Let the initial value of b be 1 and select the initial value of a such that the square of the number is as close to N as possible so as to get the least absolute value of k. Note that k may be either positive or negative but not zero.
2. Now select an integer m such that k divides a + bm and the absolute value of m2 – N is minimum.
3. Now replace a by the absolute value of (am + Nb)/k and b by the absolute value of (a + bm)/k.
4. Now replace k by a2 – Nb2. Note that k may be either positive or negative but not zero.
5. If the new value of k is 1 we have got the solution. The value of x is a and that of y is b.
6. If the value of k is not equal to 1 we go for the first iteration and repeat steps 1 to 5. We keep on doing this until the value of k becomes 1 and we thus get the values of x and y which satisfies “Pell’s equation” for the given value of N.
7. Using Brahmagupta’s lemma, if the value of k assumes a value of -1 or ± 2 for any value of a and b or if the value of k assumes a value of ± 4 and either a or b is an even number, we can stop further iterations and compose the triple with itself to get the final solution.
a2 – Nb2 = k ==> {(a2 + Nb2)/k}2 – N{(ab)/k}2 = 1

An example with N = 53 is given below.

1. Since N = 53, we initially select a = 7 and b = 1 to get k = -4.
2. We have to now select an integer m such that 4 divides 7 + m and the absolute value of m2 – 53 is minimum. Therefore m should be of the type 4t + 1 for integer values of t. To minimise the absolute value of m2 – 53, m should be 5.
3. Hence, in the first iteration we get the new values of a and b as the absolute value of (am + Nb)/k and (a + bm)/k respectively. Thus we get a = (7×5 + 53×1)/4 = 22, b = (7 + 1×5)/4 = 3 and k = (52 – 53)/(-4) = 7.
4. In the second iteration we have to select an integer m such that 7 divides 22 + 3m and the absolute value of m2 – 53 is minimum. Therefore m should be of the type 7t + 2 for integer values of t. To minimise the absolute value of m2 – 53, m should be 9.
5. The new values of a and b are now the absolute value of (am + Nb)/k and (a + bm)/k respectively. Thus we get a = (22×9 + 53×3)/7 = 51, b = (22 + 3×9)/7 = 7 and k = (92 – 53)/7 = 4.
6. In the third iteration we select m = 7 to get the revised values of a = 182, b = 25 and k = -1.
7. Since the value of k is now -1, we can use Bhaskara’s lemma an compose the triple with itself. Since a = 182, b = 25 and k = -1, we get x = (1822 + 53×252)/1 = 66249, y = 2x182x25/1 = 9100 and k = 1.
8. Thus, the solution to the equation is 662492 – 53 x 91002 = 1.
9. From this we can also get the value of √53 as 7.28010989010989. This value differs from the value obtained from the square root function of the spreadsheet by just 0.000000011 %.
10. This method always terminates with a solution (when the value of k becomes 1) for all values of N (other than perfect squares).
Project Euler 66 Solution last updated

## Discussion

### 6 Responses to “Project Euler 66 Solution”

1. Just one (dumb) quick question.
Why D has to be prime?

Posted by Tony | March 18, 2016, 7:23 AM
2. Will you please explain , how you implement these steps ??

p = k * (p/k+1) – p
p = p – int((p – sd)/k) * k

Posted by Mayank | May 29, 2013, 5:30 AM
• I added a discussion of the algorithm written by Gopal that explains all the steps required.

Posted by Mike | June 6, 2014, 8:00 PM
• Hi Mike,

Sure, it is obvious that these two serve the purpose of ensuring that a + bm (or x1 + py, in your notation) is divisable by k and that m*m – N is minimal. However, it’s really interesting how you made it happen without referring to x1 and y.

If you could explain that, I’d very much appreciate it.

Thanks,
Janko

Posted by Janko | May 27, 2015, 2:46 AM
• I solved this problem over a decade ago and then discovered the “Chaka Lacka Boom Boom” (chakravala) method to help improve performance.

The two statements you referenced are a direct result of implementing this algorithm and the following reference provides a good explanation of how it works:

https://en.wikipedia.org/wiki/Chakravala_method

Thanks

Posted by Mike Molony | June 3, 2015, 10:42 PM