We’re solving x^2 – Dy^2 = N. We’ve gotten the very easy cases out of the way: N = 0, D = 0. We’ve gotten the easy cases out of the way: D < 0, D = d^2. This leaves only the situation where D is a non-square positive integer and N is nonzero. As we will see, this one is legitimately complicated. A Pell equation is one where N = 1. Such an equation always has a solution. Now let’s suppose (x,y) is a solution to x^2 – Dy^2 = N and (r,s) is a solution to x^2 – Dy^2 = 1. If u + v sqrt(D) = (r + s sqrt(D))(x + y sqrt(D)), (u,v) will also be a solution to x^2 – Dy^2 = N. In other words, you can combine a solution of the original equation with a solution to the N = 1 equation to create a new solution of the original equation. And you can do this forever: if a solution exists, there is an infinite number of them. Luckily, it turns out that those infinite solutions are organized: there is a finite number of fundamental solutions, and each of those can generate an infinity of solutions by combining with the N = 1 solution. Our task, then, is to (a) find the N = 1 solution (the N = -1 solution is important as well, if it exists), (b) find the fundamental solutions, and (c) combine them. We can do (c) very easily using the formula above; the challenges are (a) and (b). But first, there’s an algorithm called the PQa algorithm, and there are various equivalent formulations of it. The one I will show is the one from John P. Robertson’s paper:
We’re given P_0, Q_0, and D such that Q_0 ≠ 0, D > 0 is not a perfect square, and (P_0)^2 = D (mod Q_0). (I’m assuming that, if you’re solving the Pell equation, you understand modular arithmetic, but this just means that (P_0)^2 % Q_0 == D % Q_0.) Starting with:
A_-2 = 0, A_-1 = 1
B_-2 = 1, B_-1 = 0
G_-2 = -P_0, G_-1 = Q_0
For each step i ≥ 0:
a_i = floor((P_i + sqrt(D))/Q_i)
A_i = a_i A_(i-1) + A_(i-2)
B_i = a_i B_(i-1) + B_(i-2)
G_i = a_i G_(i-1) + G_(i-2)
For each step i ≥ 1:
P_i = a_(i-1) Q_(i-1) – P_(i-1)
Q_i = (D – (P_i)^2)/Q_(i-1)
This is the skeleton of the PQa algorithm. You just update some variables in each step. I’m not going to prove anything about it. Now, just to be clear, this is how solution composition works: if (x,y) • (r,s) = (u,v), then u + v sqrt(D) = (x + y sqrt(D))(r + s sqrt(D)). Simplifying:
(x,y) • (r,s) = (u,v) = (xr + Dys,xs + yr)
If (x,y) solves the N = a equation and (r,s) solves the N = b equation, (u,v) solves the N = ab equation. This can be seen algebraically pretty easily, so I won’t go into it and return to the solution of the Pell-type equation at hand.
Since the N = ±1 case is critical, we do this first. To solve x^2 – Dy^2 = 1, we run the PQa algorithm with P_0 = 0 and Q_0 = 1, and stop when Q_k = 1 again. If k is even, then (G_(k-1),B_(k-1)) is the generator solution to N = 1, but N = -1 has no solutions. If k is odd, then (G_(k-1),B_(k-1)) is the fundamental solution to N = -1. Compose it with itself to get the generator solution to N = 1. Note that (1,0) is always the fundamental solution for N = 1, but it’s trivial. The solution you find with the algorithm is the generator solution, which you can compose with other solutions to yield new ones.
For other N, the situation is more complicated. If N = n^2, then (n,0) is a trivial solution, but there are usually others. There are multiple methods of doing this; the one I’m using is the LMM algorithm. First, find all integers f > 0 such that m = N/f^2 is an integer (it can be negative if N is negative, of course). This can be done by factoring (as in Part 1). Take the prime factorization, and replace each prime’s exponent a with floor(a/2). This is the largest square factor of N. The numbers f are the factors of this number. Now, for each m, find all z such that -|m|/2 < z ≤ |m|/2 and z^2 = D (mod |m|). This can be done by brute force, since algorithms for taking square roots in modular arithmetic are rather finicky, usually. For each z, then, run PQa with P_0 = z, Q_0 = |m|, and continue until either Q_k = ±1 (with k > 0) or one period has gone by.
How do you know when one period has gone by? A period begins at the first i where (P_i + sqrt(D))/Q_i is “reduced”, which in this context means that (P_i + sqrt(D))/Q_i > 1 and -1 < (P_i – sqrt(D))/Q_i < 0. You know the period has gone by when P_j = P_i and Q_j = Q_i. If this happens, there are no solutions for this m and z.
If Q_k = ±1, then look at r = G_(k-1), s = B_(k-1). If r^2 – Ds^2 = m, we have a set of solutions (±fr,±fs). If r^2 – Ds^2 = -m, compose it with the N = -1 solution if that exists. If it doesn’t, there are no solutions for this m and z.
There’s one problem here: we want x and y to be nonnegative, so what do we do about those ± signs in (±fr,±fs)? Two solutions are said to be equivalent if one can be generated from the other. Changing both signs will yield an obviously equivalent solution, since (-1,0) is a solution to the N = 1 equation (though not one we’re interested in since x is negative). Therefore, we can consider (±fr,fs). (fr,fs) will be its own fundamental, but (-fr,fs) can be composed with the generator solution to N = 1 until both x and y are the same sign.
We’ve ended up with a list of solutions, but the way I explained it, there might be redundant solutions. My way around this is to check for equivalence and eliminate the smaller solution when two are equivalent. (x,y) and (r,s) are equivalent if (xr – Dys)/N and (xs – yr)/N are both integers. I think you can also stop the PQa algorithm at step 0 if the end condition has been reached and get the trivial solution, if there is one, but I’m not sure there won’t be other redundant solutions so you might as well check.
Anyway, hey, we’ve solved the general Pell-type equation! It shouldn’t be too difficult to implement this procedure.
As an example, let’s work through x^2 – 5y^2 = 4.
First, we solve the N = ±1 case. This is PQa with P_0 = 0 and Q_0 = 1.
-2: A = 0, B = 1, G = 0
-1: A = 1, B = 0, G = 1
0: P = 0, Q = 1, a = 2, A = 2, B = 1, G = 2 (Q = 1, but since this is step 0, the solution it refers to is (1,0), the trivial solution)
1: P = 2, Q = 1 — this means (2,1) is a solution to N = -1, and indeed 2^2 – 5*1^2 = -1. If we compose it with itself, we get (9,4) as the generator solution for N = 1.
So we move on to the N = 4 case! We can have f = 1, m = 4 or f = 2, m = 1. For f = 1, we want to find z such that z^2 = 5 (mod 4). 5 = 1 (mod 4), so z = 1 and z = -1 both work. Let’s start with z = 1:
-2: A = 0, B = 1, G = -1
-1: A = 1, B = 0, G = 4
0: P = 1, Q = 4, a = 0, A = 0, B = 1, G = -1
1: P = -1 Q = 1 — this means (-1,1) is a solution to x^2 – 5y^2 = -4. We want +4, so we compose (-1,1) as well as (1,1) with (2,1), the solution to N = -1. In the first case, we get (3,1) as a solution, and in the second case, (7,3). Now, z = -1:
-2: A = 0, B = 1, G = 1
-1: A = 1, B = 0, G = 4
0: P = -1, Q = 4, a = 0, A = 0, B = 1, G = 1
1: P = 1, Q = 1 — this means (1,1) is a solution, but we’ve already worked with this (a computer probably wouldn’t catch this). We move on to f = 2, m = 1. We want to find z such that z^2 = 5 (mod 1). The only number mod 1 is 0, since everything is divisible by 1, so z = 0 is what we use. However, this is completely equivalent to the N = 1 case (if we were a computer, we wouldn’t know that either, by the way), and that admits the trivial solution (1,0), so this gives us the solution (2,0) (because we scale by f = 2). In a real computer program we’d want to check whether any solutions are equivalent, which we did here implicitly.
The solutions to x^2 – 5y^2 = 4 can be completely described by the fundamental solutions (2,0), (3,1), and (7,3), as well as the generator (9,4), which is not a solution (it solves N = 1, not N = 4) but can be combined any number of times with each of the fundamental solutions to yield new ones.
There we go! Wasn’t that hard, right? One thing to watch out for is that the solutions can get very, very large. As the paper shows, the generator solution for D = 991 is (379516400906811930638014896080, 12055735790331359447442538767), while the generator solution for D = 992 is a somewhat tamer (63,2). You can imagine (well, you can calculate) how large the second solution for D = 991, N = 1 is, since it’s that enormous generator composed with itself. So if you want to solve arbitrary Pell-type equations, make sure you can handle large numbers.
(Return to Part 1)