## 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.```
from Euler import sqrt, prime_sieve
def pell(d):
p, k, x1, y, sd = 1, 1, 1, 0, sqrt(d)
while k != 1 or y == 0:
p = k * (p/k+1) - p
p = p - int((p - sd)/k) * k
x = (p*x1 + d*y) / abs(k)
y = (p*y + x1) / abs(k)
k = (p*p - d) / k
x1 = x
return x
L = 1000
print "Project Euler 66 Solution:", max((pell(d),d) for d in prime_sieve(L))
```

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

Slowly swipe from either end beginning with the white vertical bar to get an idea of the starting or ending digits. For less drama, just double click the answer area. The distance between the two bars will give you an idea of the magnitude. Touch devices can tap and hold the center of the box between the two bars and choose*define*to reveal the answer.

#### 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