Solving the generalized Pell equation
Continuing the work from last week, I was aiming to solve quadratic Diophantine equations with delta = B**2 – 4*A*C > 0. I found two better references on this case which can be found in  and . These two papers describe all the necessary algorithms for solving the generalized Pell equation, to which the quadratic Diophantine equation reduces in the case delta > 0. So let’s take a look at the generalized Pell equation and methods I used to solve the equation.
The generalized Pell equation is an equation of the form x**2 – D * y**2 = N. The main case of the equation is D > 0, D is not a perfect square and abs(N) > 0. Other cases are fairly easy and solutions are straightforward.
Case D < 0
Since D < 0, equation becomes x**2 + |D| * y**2 = N, which is always positive or zero. So there are no solutions when N < 0 and only the trivial solution x = 0 and y = 0 when N = 0. If N > 0, we can note that |D| * y**2 <= N which implies y <= sqrt(N/|D|), so a brute force search on this interval would find us all the solutions. That’s what I did in this case. But there are more better ways like using methods of binary quadratic forms for positive definite forms. I hope to implement it in the near future.
Case D = 0
Equation reduces to x**2 = N, solution procedure is straight forward and no solutions exist if N < 0 or N is not a perfect square.
Case D > 0
If D is a square, then let D = r**2 for some integer r. Then the LHS is x**2 – r**2 * y**2 = N –> (x – r*y) (x + r*y) = N. Now if N = 0 either x – r*y = 0 or x+ r*y = 0 and the parametric solutions for the above two cases are (r*t, t) and (-r*t, t) respectively. If N does not equal zero, then both x – r*y and x + r*y should be factors of N. We should find an upper bound M for y so that we can search through the solutions in the range 0 <= y <= M (We don’t have to consider negative solutions, reason is explained later in the article). Let x – r*y = a and x + r*y = b. Solving for y we get y = (b – a)/(2*r). Now N = (x – r*y) * (x + r*y) = a * b so b <= N and a >= 1 since a is a positive factor o f N. which implies b <= N and -a <= -1. Now adding up the two inequalities gives b – a <= N – 1 which gives y <= (N – 1)/ (2*r) which is the required upper bound.
Now suppose D is not a perfect square. If N = 0, x**2 – D*y**2 = 0 –> x**2 = D * y**2. Suppose y is non zero. Then y should divide x (Otherwise y**2 can’t divide x**2 which is the case here). Putting x = k*y on the LHS one can see that D = k**2 which is in contradiction to the fact that D is not a perfect square. So y = 0 and consequently x = 0. Hence if N = 0 we will have only the trivial solution x = 0 and y = 0. Now if N is non zero comes our main case.
This algorithm computes the continued fraction expansion of the quadratic irrational (P + sqrt(D))/Q and some auxiliary variables, i.e of numbers like (2 + sqrt(3))/5 and (-34 + sqrt(19))/24. Note that D should not be perfect square otherwise the expression (P + sqrt(D))/Q will not be a quadratic irrational. For more details on algorithm refer . Algorithm takes three variables P, Q, D as input and output six quantities P_i, Q_i, a_i, A_i, B_i and G_i for each time it’s called for a given set of initial values (This is implemented as a generator, please see python documentation). a_i is the ith term of continued fraction expansion of (P + sqrt(D))/Q and others are required in the LMM algorithm described below. Here i >= 0.
It is also important to know that continued fraction representation of a quadratic irrationality is periodic (most of the time) after a few terms. In the solution we need to know the sum of the lengths of the aperiodic part and the periodic part. I implemented a separate function, length() to find this.
LMM algorithm finds primitive solutions of the each equivalence class of the solution set of x**2 – D * y**2 = N. We have to construct the general solution of each class afterwards. Here is a rough description how the algorithm works.
1. Make a list of f > 0 such that f**2 divides N
2. Set m = N / f**2 and find all z such that z**2 ≡ D( mod |m|) and |m|/2 < z <= |m|/2
3. Apply PQa algorithm with P = z, Q = |m| and D = D
4. Continue till there is an i >= 1 such that |Q_i| = 1 or reach the end of first period of continued fraction expansion without such a i, i.e continue till i exceeds the value returned by length()
5. If you reached an i with |Q_i| = 1 then set r = G_i-1 and s = B_i-1. if r**2 – D * s**2 = m, then add x = f*r and y = f*s to the list of solutions. Otherwise r**2 – D * s**2 = -m. Now if x**2 – D * y**2 = -1 has solutions, find it’s minimal positive solution (u, v) and add x = f*(r*u+ s*v*d) and y = f*(r*v + s*u) to the list of solutions. If x**2 – D * y**2 = -1 does not have solution, continue with next z
6.Repeat this for every z for corresponding f, and then for every f corresponding to N
We can solve the x**2 – D * y**2 = -1 by reapplying LMM algorithm but in the Diophantine module I separately solved the equation for |N| = 1 and used the solutions returned by them in the LMM algorithm.
If you are interested, a lot of cool information about solving process and the history of the Pell equation can be found in .
 Solving the equation ax^2 + bxy + cy^2 + dx + ey + f = 0, John P.Robertson, May 8, 2003. http://www.jpr2718.org/ax2p.pdf
 Solving the generalized Pell equation x**2 – D*y**2 = N, John P. Robertson, July 31, 2004. http://www.jpr2718.org/pell.pdf
 H.W Lenstra, Jr. http://www.ams.org/notices/200202/fea-lenstra.pdf