This note gives algorithms to find integer solutions x, y, to Pell equations, x2 - dy2 = N, for d a positive integer, not a square, and N a nonzero integer. What is presented here is sufficient for cases where N is 'small.' To solve other cases, you may need efficient algorithms to solve the equation x2 º d (mod |m|) where |m| is large, and to factor large integers. We give references for such algorithms, but not the algorithms themselves. Also, no proofs are given here, but references to proofs are given.
If you just want to solve a particular equation, download Keith Matthews' CALC, and use the function patz(d, N). If you want some algorithms for solving these equations, this is the place. If you want the theory behind these algorithms, see the references.
Methods are given here for solving x2 - dy2 = ±1, x2 - dy2 = ±4, and x2 - dy2 = N when N2 < d. For the general Pell equation (arbitrary N ¹ 0) there are five good methods:
Section headings are
Examples and references are given at the end.
This algorithm is at the heart of all of the algorithms to solve Pell equations presented here.
The input to the algorithm is three integers, d, P0, Q0, where d > 0 is not a square, Q0 > 0, and P02 º d (mod Q0). Recursively compute, for i ³ 0
ai = int(Pi+Öd)/Qi,
Pi+1 = aiQi - Pi, and
Qi+1 = (d-Pi+12)/Qi.
Also compute Gi and Bi as follows. Begin with G-2 = -P0, G-1 = Q0, B-2 = 1, and B-1 = 0. Then for i ³ 0, set Gi = aiGi-1 + Gi-2, and set Bi = aiBi-1 + Bi-2.
Sometimes one also computes Ai as A-2 = 0, A-1 = 1, and Ai = aiAi-1+ Ai-2 for i ³ 0. Then Gi = Q0Ai - P0Bi.
Note that Gi2 - dBi2 = (-1)i+1 Qi+1 Q0. This relation will be important to us because all of the methods of solution we discuss will involve setting Q0 = |N|, and finding those i so that (-1)i+1 Qi+1 = N/|N|. Then Gi, Bi will be a solution to the equation being considered. From a computational viewpoint, also note that, in some sense, Gi and Bi will typically be large, while Q0 and Qi+1 will be small. So this equation sometimes allows accurate computation of the left-hand-side when numbers on the left-hand-side exceed the machine accuracy available.
Exactly how far to carry these computations is discussed with each use below.
The sequence ai is the simple continued fraction expansion of (P0 + Öd)/Q0, and the Ai/Bi are the convergents to this continued fraction. Each of the sequences Pi, Qi, and ai is periodic from some point, although not necessarily the same point for all three.
References - NZM, Mollin, Rockett and Szusz, Cohen
To solve the equation x2 - dy2 = ±1, apply the PQa algorithm with P0 = 0 and Q0 = 1. There will be a smallest i with ai = 2a0, which will also be the smallest i > 0 so that Qi = 1. There are two cases to consider: this i is odd, or this i is even.
If this i is odd, the equation x2 - dy2 = -1 has solutions. The minimal positive solution is given by x = Gi-1, y = Bi-1. For any positive integer k, if k is odd then x = Gki-1, y = Bki-1 is a solution to the equation x2 - dy2 = -1, and all solutions to this equation with x and y positive are generated this way. If k is an even positive integer, then x = Gki-1, y = Bki-1 is a solution to the equation x2 - dy2 = 1, and all solutions to this equation with x and y positive are generated this way. The minimal positive solution to x2 - dy2 = 1 is x = G2i-1, y = B2i-1.
If the smallest i so that ai = 2a0 is even, then the equation x2 - dy2 = -1 does not have any solutions. For any positive integer k, x = Gki-1, y = Bki-1 is a solution to the equation x2 - dy2 = 1, and all solutions to this equation with x and y positive are generated this way. In particular, the minimal positive solution to x2 - dy2 = 1 is x = Gi-1, y = Bi-1.
The sequences Pj and aj are periodic with period i after the zero-th term, i.e. the first period is P1 to Pi for the sequence Pj, and a1 to ai for the sequence aj. The sequence Qj is periodic starting at the zero-th term, i.e. the first period is Q0 to Qi-1.
There are several methods to generate all solutions to either of the equations x2 - dy2 = ±1 once the minimal positive solution is known.
Consider first the equation x2 - dy2 = 1. If t, u is the minimal positive solution to this equation, then for the n-th positive solution xn + ynÖd = (t + uÖd)n and xn - ynÖd = (t - uÖd)n. While each positive solution corresponds to a positive n, these equations also make sense for n £ 0. There is a recursion xn+1 = txn + uynd, yn+1 = tyn + uxn. Another pair of recursions is (set x0 = 1, y0 = 0) xn+1 = 2txn - xn-1, yn+1 = 2tyn - yn-1. The comments in this paragraph apply whether or not the equation x2 - dy2 = -1 has solutions.
Now suppose the equation x2 - dy2 = -1 has solutions, let t, u be the minimal positive solution, and define xn, yn by the equation xn + ynÖd = (t + uÖd)n. Then also xn - ynÖd = (t - uÖd)n. If n is odd, xn, yn is a solution to the equation x2 - dy2 = -1, and if n is even then xn, yn is a solution to the equation x2 - dy2 = 1. All positive solutions to these two equations are so generated. The recursion xn+1 = txn + uynd, yn+1 = tyn + uxn also alternately generates solutions to the +1 and -1 equations. Another recursion is (set x0 = 1, y0 = 0) xn+1 = 2txn + xn-1, yn+1 = 2tyn + yn-1.
All solutions are given by taking the four choices of sign, ±xn, ±yn.
Perhaps the most succinct way to summarize the set of solutions is as follows. Let x, y be any solution to x2-dy2 = ±1. Let t, u be the minimal positive solution of x2 - dy2 = ±1. Then for some sign, ±1, and some integer n, x + yÖd = ±(t + uÖd)n. Note also that (t + uÖd)-1 = t - uÖd.
References: NZM, Mollin, Olds, Rockett and Szusz, Leveque, Rose, and many other sources not listed here. Many introductory books on number theory cover the Pell ±1 equation.
The material in this section is not used in any subsequent section.
In some ways, solutions to the equation x2 - dy2 = ±4 are more fundamental than solutions to the equation x2 - dy2 = ±1. The most interesting case is when d º 1 (mod 4), so we cover that first.
When d º 1 (mod 4), apply the PQa algorithm with d = d, P0 = 1, and Q0 = 2. There will be a smallest i > 0 so that ai = 2a0 - 1. This will also be the smallest i > 0 so that Qi = 2. The minimal positive solution to x2 - dy2 = ±4 is then x = Gi-1, y = Bi-1. If i is odd, it will be a solution to the -4 equation, while if i is even it will be a solution to the +4 equation and the -4 equation will not have solutions.
Periodicity of the sequences Pi, Qi, and ai is similar to that for the ±1 equation.
If d º 0 (mod 4), then for any solution to x2 - dy2 = ±4, x must be even. Set X = x/2, set Y = y, and solve X2 - (d/4)Y2 = ±1. If X, Y is the minimal positive solution to this equation, then x = 2X, y = Y is the minimal positive solution to x2 - dy2 = ±4. Alternatively, one can apply the PQa algorithm with P0 = 0 and Q0 = 2. If i is the smallest index so that ai = 2a0, then the minimal positive solution is Gi-1, Bi-1.
If d º 2 or 3 (mod 4), then by considerations modulo 4 one can see that both x and y must be even. Set X = x/2, set Y = y/2, and solve X2 - dY2 = ±1. If X, Y is the minimal positive solution to this equation, then x = 2X, y = 2Y is the minimal positive solution to x2 - dy2 = ±4. Alternatively, use the PQa algorithm with P0 = 0 and Q0 = 1, but set G-2 = 0, G-1 = 2, B-2 = 2, and B-1 = 0. If i is the smallest index so that ai = 2a0, then the minimal positive solution is Gi-1, Bi-1.
As with the ±1 equation, all solutions can be generated from the minimal positive solution. Consider first the equation x2 - dy2 = 4. If t, u is the minimal positive solution to this equation, then for the n-th solution xn + ynÖd = [(t + uÖd)n]/(2n-1) and xn - ynÖd = [(t - uÖd)n]/(2n-1). We also have the recursion xn+1 = (1/2)(txn + uynd), yn+1 = (1/2)(tyn + uxn). Another recursion is (set x0 = 2, y0 = 0) xn+1 = txn - xn-1, yn+1 = tyn - yn-1.
Now suppose the equation x2 - dy2 = -4 has solutions, let t, u be the minimal positive solution, and define xn, yn by the equation xn + ynÖd = [(t + uÖd)n]/(2n-1). Then if n is odd, xn, yn is a solution to the equation x2 - dy2 = -4, and if n is even then xn, yn is a solution to the equation x2 - dy2 = 4. All positive solutions to these two equations are so generated. The recursion xn+1 = (1/2)(txn + uynd), yn+1 = (1/2)(tyn + uxn) also alternately generates solutions to the +4 and -4 equations. Another recursion is (set x0 = 2, y0 = 0) xn+1 = txn + xn-1, yn+1 = tyn + yn-1.
The set of solutions can be summarized as follows. Let t, u be the minimal positive solution of x2 - dy2 = ±4. Then for any solution to x2 - dy2 = ±4, there is a sign, ±1, and an integer n so that (x + yÖd)/2 = (±1)[(t + uÖd)/2]n.
In some ways, the equation x2 - dy2 = ±4 is more fundamental than the equation x2 - dy2 = ±1. The numbers 1 and 4 are the only N's so that, for any d, if you know the minimal positive solution to the equation x2 - dy2 = ±N, you can generate all solutions, and you can do this without solving any other Pell equation. Also, if you know the minimal positive solution to x2 - dy2 = ±4, you can generate all the solutions to x2 - dy2 = ±1. But the converse does not hold. The best that can be said as a converse is that for d not 5 or 12, the solutions to the equation x2 - dy2 = ±4 can be derived from the intermediate steps when the PQa algorithm is used to solve the equation x2 - dy2 = ±1.
When d º 1 (mod 4), considerations modulo 4 show that for any solution to x2 - dy2 = ±4, x and y are both odd or both even. If the minimal positive solution has both x and y even, then all solutions have both x and y even. In this case, every solution to x2 - dy2 = ±1 is just one-half of a solution to x2 - dy2 = ±4. If the minimal positive solution to x2 - dy2 = ±4 has both x and y odd, then d º 5 (mod 8), every third solution has x and y even, and all other solutions have x and y odd. In this case, every solution to x2 - dy2 = ±1 is just one-half of one of the solutions to x2 - dy2 = ±4 that has both x and y even. When d º 1 (mod 4), the equation x2 - dy2 = -4 has solutions if and only if the equation x2 - dy2 = -1 has solutions.
When d º 0 (mod 4), considerations modulo 4 show that for any solution to x2 - dy2 = ±4, x is even. If the minimal positive solution has y even, then all solutions have y even (and x is always even). In this case, every solution to x2 - dy2 = ±1 is just one-half of a solution to x2 - dy2 = ±4. If the minimal positive solution to x2 - dy2 = ±4 has y odd, then every other solution has y even, and every other solution has y odd. In this case, every solution to x2 - dy2 = ±1 is just one-half of one of the solutions to x2 - dy2 = ±4 that has x and y both even. When d º 0 (mod 4), it is possible for there to be solutions to x2 - dy2 = -4, but not solutions to x2 - dy2 = -1. This happens for d = 8, 20, 40, 52 and many more values. Of course, x2 - dy2 = -1 never has solutions when d º 0 (mod 4).
When d º 2 (mod 4) or d º 3 (mod 4), all solutions to x2 - dy2 = ±4 have both x and y even. Every solution to x2 - dy2 = ±1 is just one-half of a solution to x2 - dy2 = ±4. The equation x2 - dy2 = -4 has solutions if and only if the equation x2 - dy2 = -1 has solutions.
References - Cohen, NZM, Mollin, Leveque. Cohen treats the cases d º 1 (mod 4) for d squarefree, and d = 4r for r º 2 or 3 (mod 4), r squarefree . The above material is not really addressed directly in either of NZM or Mollin. But the only matter above that is not trivially derived from material in one or both of these sources is the proof that the method for solving the equation works in the case d º 1 (mod 4). Here, one can imitate the proof in NZM for the equation x2 - dy2 = ±1, and make use of Mollin's theorem 5.3.4, p. 246. This will result in a proof for all d º 1 (mod 4), d not a square; not just for the d treated in Cohen. Leveque only treats the generation of all solutions from the base solution.
If r, s is a solution to x2 - dy2 = N, and t, u is any solution to x2 - dy2 = 1, then for x = rt + sud, y = ru + st, x, y is a solution to x2 - dy2 = N. This follows from the relation (r2 - ds2)(t2 - du2) = (rt + sud)2 - d(ru + st)2. This fact can be used to separate solutions to x2 - dy2 = N into equivalence classes. Two solutions x, y and r, s are equivalent if there is a solution t, u to t2 - du2 = 1 so that x = rt + sud and y = ru + st.
It may help to view the set of solutions geometrically. If N > 0, then, as an equation in real numbers, x2 - dy2 = N is a hyperbola with the x-axis as its axis, and the y-axis as an axis of symmetry. The asymptotes are the lines x ± yÖd = 0. Let t, u be the minimal positive solution to x2 - dy2 = 1. Draw the graph of x2 - dy2 = N over the reals. Mark the point (ÖN, 0), which is on this graph. Now mark the point (tÖN, uÖN), which is also on the graph. Continue marking points so that if (x, y) is the most recent point marked, then the next point marked is (xt + yud, xu + yt). All of the points marked so far, apart from the first, have x > 0 and y > 0. Now, for each point (x, y) that has been marked, mark all of the points (±x, ±y) not yet marked.
The marked points divide the graph into intervals. Make the interval ((ÖN, 0), ( tÖN, uÖN)] a half-open interval, and then make the other intervals on this branch half-open by assigning endpoints to one interval. Make the intervals on the other branch half-open by mapping (x, y) in the right branch to (-x, -y) on the left branch. If there are integer solutions to x2 - dy2 = N, then
1) No two solutions within the same (half-open) interval are equivalent,
2) Every interval has exactly one solution in each class, and
3) The order of solutions by class is the same in every interval.
Instead of starting with the point (ÖN, 0), we could have started with any point (r, s) on the graph, and marked off the points corresponding to (r + sÖd)*(±1)*(t + uÖd)n. The above three comments would still apply.
The situation is similar in N < 0, except that the graph has the y-axis as its axis, and the x-axis is an axis of symmetry.
If x2 - dy2 = -1 has solutions, then any of these solutions can be used to form a correspondence between solutions to x2 - dy2 = N and -N.
Within a class there is a unique solution with x and y nonnegative, but smaller than any other nonnegative solution. This is the minimal nonnegative solution for the class. There is also either one or two solutions so that y is nonnegative, and is less than or equal to any other nonnegative y in any solution x, y within the class. If there is one such solution, it is called the fundamental solution. If there are two such solutions, then they will be equivalent and their x-values will be negatives of each other. In this case, the solution with the positive x-value is called the fundamental solution for the class. For N > 0, the fundamental solutions are on the hyperbola in the intervals
(ÖN, 0) to (Ö((r + 1)N/2), Ö(N(r - 1)/(2d))), and
(-ÖN, 0) to (-Ö((r + 1)N/2), Ö(N(r - 1)/(2d))).
For the first interval, the endpoints should be included, while for the second interval they should be excluded.
For N < 0, the fundamental solutions are in the interval
(-Ö((r - 1)|N|/2), Ö(|N|(r + 1)/(2d))) to
(Ö((r - 1)|N|/2), Ö(|N|(r + 1)/(2d))),
with midpoint (0, Ö(-N/d)). In this interval, the first point should be excluded, and the last point included.
When tabulating solutions, it is usually convenient to make a list consisting of one solution from each class. Often, this list will be either the minimal nonnegative solutions, or the fundamental solutions. Given any solution in a class, it is easy to find the fundamental solution or the minimal nonnegative solution for that class.
To summarize, given any solution in a class, all solutions in that class are found by applying solutions to the equation x2 - dy2 = 1. If r, s is any particular solution to x2 - dy2 = N, x, y is any other solution to the same equation in the same class as r, s, and if t, u is the minimal positive solution to the equation x2 - dy2 = 1, then for some choice of sign, ±1, and for some integer n, x + yÖd = ±(r + sÖd)(t + uÖd)n.
There are recursion relations among solutions similar to those presented for the ±1 and ±4 equations. For instance, if (x1, y1), (x2, y2), and (x3, y3) are three solutions in the same class, in consecutive intervals, and t, u is the minimal positive solution to x2 - dy2 = 1, then x3 = 2tx2 - x1 and y3 = 2ty2 - y1.
References - NZM, Mollin, Chrystal, Leveque, Rose
When 1 < N2 < d, apply the PQa algorithm with d = d, P0 = 0, Q0 = 1. Continue the computations until you reach the first i > 0 with Gi2 - dBi2 = 1 (i.e. Qi+1 = 1 and i+1 is even). For 1 £ j £ i, if Gj2 - dBj2 = N/f2 for some f > 0, add fGj, fBj to the list of solutions. When done, the list of solutions will have the minimal positive member of each class.
The list of all solutions can be generated using the methods of the previous section. Alternatively, all positive solutions can be generated by extending the PQa algorithm indefinitely.
References - NZM, Mollin, Chrystal
Let t, u be the minimal positive solution to x2 - dy2 = 1. If N > 0, set ylim1 = 0, and ylim2 = Ö (N(t - 1)/(2d)). If N < 0, set ylim1 = Ö((-N)/d), and ylim2 = Ö((-N)(t + 1)/(2d)). For ylim1 £ y £ ylim2, if N+dy2 is a square, set x = Ö(N + dy2). If (x, y) is not equivalent to (-x, y), add both to the list of solutions, otherwise just add (x, y) to the list. When finished, this list gives the fundamental solutions.
This method works well if ylim2 is not too large, which means that Ö(|N|(t ± 1)/(2d)) is not too large. You must be able to perform the search between the limits ylim1 and ylim2.
To generate all solutions from these, see the section Structure of solutions to x2 - dy2 = N .
References - Mollin, Leveque, Rose
Make a list of f > 0 so that f2 divides N. For each f in this list, set m = N/f2. Find all z so that -|m|/2 < z £ |m|/2 and z2 º d (mod |m|). For each such z, apply the PQa algorithm with P0 = z, Q0 = |m|, d = d. Continue until either there is an i with Qi = ±1, or, without having reached an i with Qi = ±1, you reach an i so that there is a j with j < i, Pi = Pj, and Qi = Qj. In the latter case, there will not be any i with Qi = ±1. If you reached an i with Qi = ±1, then look at r = Gi-1, s = Bi-1. If r2 - ds2 = m, then add x = fr, y = fs to the list of solutions. Otherwise, r2 - ds2 = -m. If the equation t2 - du2 = -1 has solutions, let the minimal positive solution be t, u, and add x = f(rt + sud), y = f(ru + st) to the list of solutions.
When you have done every f, and every z for each f, the list of solutions will have one member from each class. These solutions will be either fundamental or the minimal positive solution for the class.
To generate all solutions from these, see the section Structure of solutions to x2 - dy2 = N . Alternatively, for each z that gives rise to solutions, you can extend the PQa algorithm indefinitely.
If |N| is large, it may be necessary to have an efficient method to factor N to make the list of f's so that f2 divides N. The literature on factoring is vast. Many mathematical software packages, such as Maple or PARI, have efficient factoring systems built in. Methods for factoring integers include Pollard's number field sieve, the special number field sieve, Pomerance's quadratic sieve, the Brillhart-Morrison continued fraction factoring algorithm, D. Shanks' square-free factorization, Pollard's rho method, and the elliptic curve method. See NZM, Pomerance, Bressoud, Mollin, and many other sources.
Also, if |N| is large, it may be necessary to have an efficient method to solve the equation x2 º d (mod |m|). Cohen gives some methods for solving x2 º d (mod p) where p is an odd prime. From such solutions, one can readily solve the more general equation x2 º d (mod |m|).
When N is large, the other methods (Lagrange's system of reductions, cyclic method, binary quadratic forms) also require efficient methods to factor integers and to solve x2 º d (mod m).
Keith Matthews has a program CALC that applies this method.
I call this the LMM method because it apparently has been independently discovered by Lagrange, Matthews, and Mollin. This presentation is based on Matthews' paper.
Reference - Matthews
Apply the PQa algorithm with d = 109, P0 = 0 and Q0 = 1. The following table gives the index, i, and then the several calculated quantities for i = -2 to 30.
i Pi Qi ai Gi Bi G2 - 109B2 -2 0 1 -109 -1 1 0 1 0 0 1 10 10 1 -9 1 10 9 2 21 2 5 2 8 5 3 73 7 -12 3 7 12 1 94 9 7 4 5 7 2 261 25 -4 5 9 4 4 1138 109 15 6 7 15 1 1399 134 -3 7 8 3 6 9532 913 3 8 10 3 6 58591 5612 -15 9 8 15 1 68123 6525 4 10 7 4 4 331083 31712 -7 11 9 7 2 730289 69949 12 12 5 12 1 1061372 101661 -5 13 7 5 3 3914405 374932 9 14 8 9 2 8890182 851525 -1 15 10 1 20 181718045 17405432 9 16 10 9 2 372326272 35662389 -5 17 8 5 3 1298696861 124392599 12 18 7 12 1 1671023133 160054988 -7 19 5 7 2 4640743127 444502575 4 20 9 4 4 20233995641 1938065288 -15 21 7 15 1 24874738768 2382567863 3 22 8 3 6 169482428249 16233472466 -3 23 10 3 6 1041769308262 99783402659 15 24 8 15 1 1211251736511 116016875125 -4 25 7 4 4 5886776254306 563850903159 7 26 9 7 2 12984804245123 1243718681443 -12 27 5 12 1 18871580499429 1807569584602 5 28 7 5 3 69599545743410 6666427435249 -9 29 8 9 2 158070671986249 15140424455100 1 30 10 1 20 3231012985468390 309474916537249 -9We have a0 = 10, and the first i so that ai = 2a0 is i = 15, at which point a15 = 20. Whence the period of ai is 15, which is odd, and so the equation x2 - 109y2 = -1 has solutions. The minimal positive solution to x2 - 109y2 = -1 is x = 8890182, y = 851525. The minimal positive solution to x2 - 109y2 = 1 is x = 158070671986249, y = 15140424455100.
x2 - 109y2 = ±4.
As d º 1 (mod 4), apply the PQa algorithm with d = 109, P0 = 1 and Q0 = 2. The following table gives the index, i, and the standard quantities for i = -2 to 14.
i Pi Qi ai Gi Bi G2 - 109B2 -2 -1 1 -108 -1 2 0 4 0 1 2 5 9 1 -28 1 9 14 1 11 1 12 2 5 6 2 31 3 -20 3 7 10 1 42 4 20 4 3 10 1 73 7 -12 5 7 6 2 188 18 28 6 5 14 1 261 25 -4 7 9 2 9 2537 243 28 8 9 14 1 2798 268 -12 9 5 6 2 8133 779 20 10 7 10 1 10931 1047 -20 11 3 10 1 19064 1826 12 12 7 6 2 49059 4699 -28 13 5 14 1 68123 6525 4 14 9 2 9 662166 63424 -28We have a0 = 5, and the first i so that ai = 2a0 - 1 is i = 7, at which point a7 = 9. Whence the period of ai is 7, which is odd, and so the equation x2 - 109y2 = -4 has solutions. The minimal positive solution to x2 - 109y2 = -4 is x = 261, y = 25. The minimal positive solution to x2 - 109y2 = 4 is x = 68123, y = 6525.
Note that the third solution to x2 - 109y2 = ±4 can be computed from [(261 + 25Ö(109))3]/4 = 17780364+1703050Ö(109). Upon dividing by 2, we get the minimal positive solution x = 8890182, y = 851525 to x2 - 109y2 = -1.
x2 - 129y2 = -5
Here N2 < d. Apply the PQa algorithm with d = 129, P0 = 0, Q0 = 1.
i Pi Qi ai Gi Bi G2 - 129B2 -2 0 1 -129 -1 1 0 1 0 0 1 11 11 1 -8 1 11 8 2 23 2 13 2 5 13 1 34 3 -5 3 8 5 3 125 11 16 4 7 16 1 159 14 -3 5 9 3 6 1079 95 16 6 9 16 1 1238 109 -5 7 7 5 3 4793 422 13 8 8 13 1 6031 531 -8 9 5 8 2 16855 1484 1 10 11 1 22 376841 33179 -8The only f > 0 so that f2 divides -5 is f = 1. Reviewing the above for Gj2 - 129Bj2 = -5, we find solutions x, y equal to 34, 3 and 1238, 109. Thus, there are two classes of solution, and these are the minimal positive solutions for these classes.
x2 - 61y2 = 15 by brute-force search
The minimal positive solution to x2 - 61 y2 = 1 is x = 1,766,319,049, y = 226,153,980. As N = 15 is positive, the lower search limit for y is 0, and the upper limit is Ö((15*(1766319049 -1))/(2*61)), which is about 14,736.702. So we search on y from 0 to 14,736. Only y = 11 and y = 917 yield solutions, so the four fundamental solutions are x = ±86, y = 11, and x = ±7162, y = 917.
x2 - 61y2 = 15 by LMM
The only f > 0 so that f2 divides 15 is f = 1. Set m = 15. The z's with -15/2 < z £ 15/2 and z2 º 61 (mod 15) are z = ±1, z = ±4.
Upon performing the PQa algorithm with P0 = 1, Q0 = 15, and d = 61, the first Qi = ±1 occurs at Q9 = 1. The corresponding solution has x = G8 = 2593 and y = B8 = 332. For this x, y, x2 - 61y2 = -15. The equation x2 - 61y2 = -1 has solutions, and the minimal positive solution is x = 29718, y = 3805. Applying this solution to the solution 2593, 332 gives the solution x = 154117634, y = 19732741 to the equation x2 - 61y2 = 15. This is equivalent to the fundamental solution -86, 11.
Performing the PQa algorithm with P0 = -1, Q0 = 15, and d = 61, gives the first Qi = ±1 at Q4 = 1, yielding the fundamental solution 86, 11 to the equation x2 - 61y2 = 15.
Performing the PQa algorithm with P0 = 4, Q0 = 15, and d = 61, gives the first Qi = ±1 at Q10 = 1, yielding the fundamental solution 7162, 917 to the equation x2 - 61y2 = 15.
Performing the PQa algorithm with P0 = -4, Q0 = 15, and d = 61, gives the first Qi = ±1 at Q3 = 1, yielding the solution 31, 4 to the equation x2 - 61y2 = -15. Applying the minimal positive solution to x2 - 61y2 = -1 gives the solution x = 1849678, y = 236827 to the equation x2 - 61y2 = 15. This is equivalent to the fundamental solution -7162, 917.
References
David M. Bressoud, Factorization and Primality Testing, Springer-Verlag, NY, 1989. Covers the Pomerance quadratic sieve factoring method, the elliptic curve factoring method, and more. Also gives the PQa algorithm for solving Pell equations x2 - dy2 = ±1.
G. Chrystal, Algebra, An Elementary Text-Book (a.k.a. Textbook of Algebra), Part II, Dover, NY, 1961. Other editions include Adam and Charles Black, 1900, and Chelsea, perhaps published in the 1950's, and AMS currently. Chapter XXXII, pages 423 to 452, begins a discussion of continued fractions. Chapter XXXIII, pp. 453 to 490, discusses recurring continued fractions. In Chapter XXXIII, Sections 15 to 17, pages 478 to 481, recurring continued fractions are studied and applied to solve the equations x2 - dy2 = ±1, and x2 - dy2 = m for m2 < d. Chapter XXXIII, Section 18, pages 482 to 485, discusses Lagrange's method of reduction for the case m2 > d, but apart from the fact that the heading on page 482 is "Lagrange's Chain of Reductions," the section does not reference Lagrange. Section 19, on page 486, discusses the cases d < 0 and d > 0, d a square. Section 20, pages 486 to 488, reduces the general equation ax2+2hxy+by2+2gx+2fy+c = 0 to a Pell equation covered previously. Data for the current AMS edition is: Algebra, an Elementary Text-Book for the Higher Classes of Secondary Schools and for Colleges: Seventh Edition - G. Chrystal - AMS | CHEL, 1964, 1212 pp., Hardcover, ISBN 0-8218-1931-3, List: $40, All AMS Members: $36, CHEL/84.H
Henri Cohen, A Course in Computational Algebraic Number Theory, Springer-Verlag, 1993. Section 1.5, pp. 31-36, covers methods to solve x2 º d (mod p), for p prime. Algorithm 5.7.2 in Section 5.7, pages 264 to 274, solves the equation x2 - dy2 = ±4 when d º 1 (mod 4) for d squarefree, or d = 4r for r º 2 or 3 (mod 4), r squarefree.
Harold M. Edwards, Fermat's Last Theorem, Springer-Verlag, NY, 1977. Chapters 7 (pp. 245 to 304) and 8 (pp. 305 to 341), especially sections 8.2 (pp. 313-318) and 8.7 (pp. 339-341) apply the cyclic method to solve equations of the form ax2 + bxy + cy2 + dx + ey + f = 0.
Adolf Hurwitz, Lectures on Number Theory, Springer-Verlag, New York, 1986. Chapter 6 is a wonderful exposition of methods to solve binary quadratic form equations, Ax2 + Bxy + Cy2 = N. Much theory of simple continued fractions is also developed.
William Judson Leveque, Topics in Number Theory, Volume 1, Addison-Wesley, New York , 1956. Chapter 8, sections 1 to 3, pages 137 to 148, discuss the equations x2 - dy2 = ±1, x2 - dy2 = ±4, and x2 - dy2 = N for general N. Limits on the size of |x| for fundamental solutions are given in Theorem 8-9, page 147, and in exercise 3, page 148 (the first inequality in this exercise should be 0 £ u, not 0 < u; limits sharper than those given in Theorem 8-9 or exercise 3 are possible).
Keith Matthews, The diophantine equation x2-Dy2 = N, D > 1, in integers, Expositiones Mathematicae, 18 (2000), 323-331. Gives the LMM method for solving x2 - dy2 = N for any nonzero N. Available at Keith Matthews' publications page.
Richard E. Mollin, Fundamental Number Theory with Applications, CRC Press, Boca Raton, 1998. Chapter 5, pages 221 to 272, discusses the continued fractions generally. In particular, Section 5.3, pages 238 to 250 studies periodic continued fractions, and applies this study to find all solutions to the Pell equation x2-dy2 = ±1 for d>0, d not a square. The PQa algorithm for computing the continued fraction expansion of a quadratic irrational is discussed in exercise 5.3.6, p. 251. While the method above for the ±4 equation is not discussed explicitly, Mollin's Theorem 5.3.4 on page 246 gives the main machinery needed to prove that that method is correct. Mollin also discusses the continued fraction expansion of (1+Öd)/2 for d º 1 (mod 4) in exercise 5.3.14 on page 252. For the general equation x2-dy2 = m, for any positive nonsquare d and any m, limits to search on x or y are given in Chapter 6 on pages 299 and 300. Theorem 5.2.5 on page 232 gives the main criterion for solving x2-dy2 = m when m2 < d. Corollary 6.2.1 (page 305) to Theorem 6.2.7 (page 302) gives the essence of Lagrange's system of reduction, which can be used to solve x2 - dy2 = m for d > 0, d not a square, m2 > d. Chrystal gives a more complete exposition of Lagrange's system of reduction.
Ivan Niven, Herbert S. Zuckerman, and Hugh L. Montgomery (NZM), An Introduction to the Theory of Numbers, Fifth Edition, John Wiley & Sons, Inc., New York, 1991. Sections 7.1 to 7.7, pages 325 to 351 cover continued fractions generally. Section 7.8, pages 351 to 356 covers the Pell ±1 equation, and the case x2 - dy2 = n where n2 < d. The PQa algorithm for computing the solution to the ±1 Pell equation is given in Section 7.9, page 358. This is also covered in Section 7.7, pages 346-348. Page 295 has notes on factoring integers.
C. D. Olds, Continued Fractions, MAA, 1963. An easy introduction to simple continued fractions, and the Pell equation x2 - dy2 = ±1. Does not develop the PQa algorithm.
Carl Pomerance, A Tale of Two Sieves, Notices of the American Mathematical Society, Vol. 43, No. 12, December 1996, pages 1473 to 1485. Discusses the Pomerance quadratic sieve factoring algorithm and the Pollard number field sieve factoring algorithm. Available as a pdf file at the AMS journals page
Andrew M. Rockett and Peter Szusz, Continued Fractions, World Scientific, 1992. Good introduction to continued fractions generally, with treatment of x2 - dy2 = ±1.
H. E. Rose, A Course in Number Theory, Clarendon Press, 1988. Chapter 7, section 3, pages 125 to 128 treats the equation x2 - dy2 = ±1. Theorem 3.3, page 128, gives limits on |x| for a fundamental solution to the general case x2 - dy2 = m (although the first inequality has to be changed from 0 < u to 0 £ u).
Last Modified May 25, 2002 John P. Robertson JPR2718@AOL.COM Copyright John P. Robertson 2000, 2001, 2002