Solving Pell's equations with the Chakravala method and Clojure 28 Jan 2012
Another day, another Project Euler problem. This time, I'm going to talk about Problem 66, which is an example of what makes Project Euler so awesome — a problem that has no brute force solution, and requires you to learn something new (and implement it in code) in order to solve.
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\times2^2 = 1\)
\(2^2 – 3\times1^2 = 1\)
\(9^2 – 5\times4^2 = 1\)
\(5^2 – 6\times2^2 = 1\)
\(8^2 – 7\times3^2 = 1\)Hence, by considering minimal solutions in \(x\) for \(D \le 7\), the largest \(x\) is obtained when \(D = 5\).
Find the value of \(D \le 1000\) in minimal solutions of \(x\) for which the largest value of \(x\) is obtained.
The subset of Diophantine equations asked about in this problem are called Pell's equations. These equations are interesting for a number of reasons, one of which being that \(\frac{x}{y}\) is a good rational approximation of \(\sqrt{n}\). Luckily for us, there are a couple of approaches that will solve the equation for any non-square \(n\), so it was just a matter of choosing one. It seems like the continued fraction approach would the fastest, but I decided instead to implement the 1000 year old Chakravala method, which was described by Bhāskara II in the 12th century (call me what you will, but implementing algorithms from the distant past on modern hardware is really fun).
First, the code (all of which is available in the 0.0.2-snapshot release of clj-numbers), then some explanation of how it works.
The Chakravala method has four steps:
- Find any integer solution for \(a^2 - nb^2 = k\) where \(k \ne 0\)
- Compose the triple \((a, b, k)\) with \((m, 1, m^2 - n)\) to get the triple \((am + nb, a + bm, k(m^2 − n))\), which is then scaled by k (according to Bhaskara's lemma).
- Choose a positive integer value for \(m\) such that \(\frac{a + bm}{k}\) is an integer, and which results in the smallest possible absolute value for \(m^2 - n\).
- Calculate a new triple using those values. If \(k = 1\), then \(a\) and \(b\) are the minimal integer solution for the equation. Otherwise, return to step 2 with the new triple.
The code above breaks each of those three steps into a function, and then runs through them until a solution is found. find-initial returns an \((a, b, k)\) triple that has the smallest possible values. find-optimal-m then returns the value for \(m\) that satisfies the requirements in step #3, and finally step takes the original \((a, b, k)\) triplet and the computed value for \(m\) and returns a new triplet. This repeats until a new triplet's \(k\) equals 1, at which point \(a\) and \(b\) are returned as the solution.
Now that we have working code to calculate the solution to Pell's equation for any non-square \(n\), the rest of the problem is just a bit of list comprehension. Here's the code for the solution:
The solution takes ~120 milliseconds to run on my computer, making me think that my implementation of the Chakravala method is reasonably efficient. One obvious optimization would be to create a more efficient implementation of find-optimal-m, which right now potentially checks a lot of numbers before it finds a solution, but that's an exercise for another day.
If you liked this, you may also be interested in my other Project Euler solutions in Clojure.
- Next: Fixing math education
- Previous: Introducing clj-numbers
