// you’re reading...
1 Star2 Stars3 Stars4 Stars5 Stars (22 votes, average: 5.00 out of 5)
Loading...

Project Euler Solutions

Project Euler 66 Solution

Project Euler 66 Solution

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.

download arrowUse 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

Post a comment