## Project Euler 66: Investigate the Diophantine equation x^{2} − Dy^{2} = 1.

#### Problem Description

Consider quadratic Diophantine equations of the form:

*x*^{2} – D*y*^{2} = 1

For example, when D=13, the minimal solution in *x* is 649^{2} – 13×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.

#### 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 a^{2} – Nb^{2} = k, can be composed with the trivial triple (m, 1, m^{2} – N) to get a new triple (am + Nb, a + bm, k(m^{2} – N)). Scaling down by k, using Bhaskara’s lemma, we get,

a^{2} – Nb^{2} = k ==> {(am + Nb)/k}^{2} – N{(a + bm)/k}^{2} = (m^{2} – N)/k

The following is the algorithm for solving “Pell’s equation”, x^{2} – Ny^{2} = 1, as given by the Indian mathematicians who were the first to solve this equation:

- 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 a
^{2}– Nb^{2}= 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. - Now select an integer m such that k divides a + bm and the absolute value of m
^{2}– N is minimum. - Now replace a by the absolute value of (am + Nb)/k and b by the absolute value of (a + bm)/k.
- Now replace k by a
^{2}– Nb^{2}. Note that k may be either positive or negative but not zero. - 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.
- 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.
- 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.

a^{2}– Nb^{2}= k ==> {(a^{2}+ Nb^{2})/k}^{2}– N{(ab)/k}^{2}= 1

An example with N = 53 is given below.

- Since N = 53, we initially select a = 7 and b = 1 to get k = -4.
- We have to now select an integer m such that 4 divides 7 + m and the absolute value of m
^{2}– 53 is minimum. Therefore m should be of the type 4t + 1 for integer values of t. To minimise the absolute value of m^{2}– 53, m should be 5. - 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 = (5
^{2}– 53)/(-4) = 7. - In the second iteration we have to select an integer m such that 7 divides 22 + 3m and the absolute value of m
^{2}– 53 is minimum. Therefore m should be of the type 7t + 2 for integer values of t. To minimise the absolute value of m^{2}– 53, m should be 9. - 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 = (9
^{2}– 53)/7 = 4. - In the third iteration we select m = 7 to get the revised values of a = 182, b = 25 and k = -1.
- 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 = (182
^{2}+ 53×25^{2})/1 = 66249, y = 2x182x25/1 = 9100 and k = 1. - Thus, the solution to the equation is 66249
^{2}– 53 x 9100^{2}= 1. - 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 %.
- 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*

Just one (dumb) quick question.

Why D has to be prime?

Stranger still is that they are all primes of the form 4k+1. A complete explanation, much better that I could hope to do, is located here:

http://mathworld.wolfram.com/PellEquation.html

Will you please explain , how you implement these steps ??

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

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

I added a discussion of the algorithm written by Gopal that explains all the steps required.

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

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