# Pell equation x² - Dy² = k

This equation is meaningfull only when searching integer solutions.

## General considerations

If D negative (x² + Ay² = k) there is a finite number of solutions, if there are.
Also if D is a square x² - A²y² = k can be written (x + Ay)(x - Ay) = K that is x + Ay and x - Ay have to be choosen among divisors of K, finitely many.

We'll study here only the case  D positive not a square  and first of all the case k = 1 :
x² - Dy² = 1 (proper Pell equation).

We then have the following results :

• Equation x² - Dy² = 1 has always the trivial solution x = 1, y = 0

• Equation x² - Dy² = 1 has always infinitely many non trivial solutions.
They can be sorted in increasing values of x or y
For two solutions x1 < x2 <=> y1 < y2
The smallest non trivial solution is named fundamental solution and noted as (x0, y0)

• If (x1, y1) and (x2, y2) are two solutions, then (x, y) = (x1, y1)o(x2, y2) given by :
x = x1x2 + Dy1y2 and y = x1y2 + x2y1 also, from Brahmagupta relation
(X² + AY²)(X'² + AY'²) = (XX' -/+ AYY')² + A(XY' +/- X'Y)²

• All solutions are deduced from the fundamental solution using previous relation repeatedly, resulting into :

xn+1 = x0xn + Dy0yn,    yn+1 = y0xn + x0yn

We can then deduce recurence relation for the sequences xn and yn :

xn+2 = 2x0xn+1-xn    yn+2 = 2x0yn+1-yn

We can also write :  xn + yn√D = (x0 + y0√ D)n+1  or the conjugate form, replacing + by -

This gives closed forms for xn and yn as a function of n. These formulaes seems to be of little use because of irrational terms in √D, although giving integral xn and yn !
(They are easy to write from the previous relation combined with its conjugate)

## Example x² - 8y² = 1

It is not too difficult here to guess a non trivial solution x = 3, y = 1. This is in fact the fundamental solution, as any lower solution would have y < 1 hence would be the trivial solution (1, 0).
Applying iteration relations, we get the solutions :
(3,1) (17,6) (99,35) (577,204) ...

## Fundamental solution

In less easy cases, remains to find the fundamental solution.
Mention the treacherous question from Fermat to Frenicle with equation x² - 61y² = 1.
Thus the fundamental solution of this equation is x0 = 1766319049 y0 = 226153980. (remind that the fundamental solution is the smallest solution other than x = 0, y = 1 ! )
The size of this solution is unrelated with the size of D, there is no simple rule between both values.
For instance the nearly similar equation x² - 63y² = 1 has an obvious fundamental solution x = 8, y = 1.
x² - 13y² = 1 needs already an efficient method for searching the fundamental solution.
Hopefully, there is one, given here as an algorithm :

 Let U0 = 0, V0 = 1, a0 = E[√D], P0 = a0, Q0 = 1 and by convention P-1 = 1 and Q-1 = 0 Noting E[x] = integer part of x. Then repeat :     Un+1 = anVn - Un,   Vn+1 = (D - Un+1²)/Vn,   an+1 = E[ (a0 + Un+1)/Vn+1 ],     (Pn+1, Qn+1) = an+1(Pn, Qn) + (Pn-1, Qn-1) until ak > a0 (or Vk = 1). Then (Pk-1, Qk-1) is solution of x² - Dy² = (-1)k. If k even it is the required solution. else, solution of x² - Dy² = 1 is obtained from solution (P, Q) of x² - Dy² = -1 by :     (P² + DQ², 2PQ)

This algorithm is a specific case of PQa algorithm giving the expansion as a continued fraction of a quadratic number (u + √D)/v with here u = 0 and v = 1.

Apply this to x² - 13y² = 1 :

 U0 = 0 V0 = 1 a0 = 3 P0 = 3 Q0 = 1 U1 = 3 V1 = (13 - 9)/1 = 4 a1 = E[(3 + 3)/4] = 1 P1 = 4 Q1 = 1 U2 = 4 - 3 = 1 V2 = (13 - 1)/4 = 3 a2 = E[(3 + 1)/3] = 1 P2 = 7 Q2 = 2 U3 = 3 - 1 = 2 V3 = (13 - 4)/3 = 3 a3 = E[(3 + 2)/3] = 1 P3 = 11 Q3 = 3 U4 = 3 - 2 = 1 V4 = (13 - 1)/1 = 4 a4 = E[(3 + 1)/4] = 1 P4 = 18 Q4 = 5 U5 = 4 - 1 = 3 V5 = (13 - 9)/4 = 1 a5 = E[(3 + 3)/1] = 6 > 3 => stop

5 being odd, (18, 5) is solution of x² - 13y² = -1.
The fundamental solution of x² - 13y² = 1 is (18² + 13×5², 2×18×5) = (649, 180) which was not obvious without this algorithm.

Let's go on with the general case k ≠ 1, and first with

## Equation x² - Dy² = -1

This equation may have no solutions, but when it has, there are infinitely many.
Here also, they can be sorted in increasing order and the smallest is named "fundamental solution" or "primitive solution" and noted as (x0, y0)
All other solutions are given by  xn + yn√D = (x0 + y0√ D)2n+1

or the conjugate form by replacing + by -

Note : solutions of equation x² - Dy² = 1 are obtained from the fundamental solution
(x0, y0) of x² - Dy² = -1 by : xn + yn√D = (x0 + y0√D)2n
The fundamental solution of x² - Dy²=1 is then (x0² + Dy0², 2x0y0).
We also have the recurrence relations :

xn+1 = (x0² + Dy0²)xn + 2Dx0y0yn,   yn+1 = 2x0y0xn + (x0² + Dy0²)yn

As well as relations in xn+2 = axn+1 + bxn that are left as an exercice.

As for the searching of the fundamental solution, just use the previous algorithm, and note that when k is even, then x² - Dy² = -1 has no solution. There is no easy criterion to know it in advance. Note the case when D is a prime p :
Equation x² - py² = -1 has solutions if and anly if p = 4k + 1.

### Example :

Equation x² - 2y² = -1
x = 1, y = 1 is an obvious solution, and is fundamental as any smaller solutions would have x and y < 1.
The fundamental solution of x² - 2y² = +1 is then x = 1² + 2*1² = 3 and y = 1² + 1² = 2 .
All other solutions are then xn+1 = 3xn + 4yn,   yn+1 = 2xn + 3yn, that is :
(1,1) (7,5) (41,29) (239,169) (1393,985) ...

## Generalized Pell equation

That is x² - Dy² = K with |K| ≠ 1
Here also, this equation may have no solutions, but if there are, there are infinitely many. On the other hand, there is no more a fundamental solution but a set of primitive solutions
Let's continue to name (x0, y0) the fundamental solution of x² - Dy² = 1.
Then all solutions of x² - Dy² = K are obtained from the set of independant primitive solutions (ξ, η) by the set of relations

xn + yn√D = (x0 + y0√D)n(ξ + η√D)  and also the corresponding recurrence relations.

The problem is now to find these independant primitive solutions. A first method is to scan systematically all values of x or y, the theory giving just an upper bound. (remind the equation may have no solutions). This is the most efficient method if the upper bound is not too big.

But at first we must explain what are independant solutions, or more precisely when they are not, that is when they are said to be "equivallent".
Two solutions (x,y) and (x',y') are equivalent if there is a solution (X,Y) of x² - Dy² = 1 with :

x + y√D = (X + Y√D)(x' + y'√D)

This criterion being not very usefull, add the following :

 Two solutions (x,y) and (x',y') of x² - Dy² = K are equivallent if and only if :    xx' - Dyy' and x'y - xy' are multiple of K.

Remains to find bounds for the primitive solutions :

 If K > 0, the primitive solutions are bounded by 0 < x ≤ √( (x0 + 1)K/2 ) or 0 ≤ y ≤ y0√ ( K/2(x0 + 1) ) If K < 0, the primitive solutions are bounded by 0 ≤ x ≤ √( (x0 - 1)|K|/2 ) or 0 < y ≤ y0√( |K|/2(x0 - 1) )

If there is no solutions in the bounds, there is no solutions at all.

We would like to find a simple criterion for x² - Dy² = K to have solutions by writing
x² ≡ K mod (D), that is K must be a quadratic residue of D.
The only thing we can say then is that when K is not a quadratic residue of D, then the Pell equation has no solutions. For instance equation x² - 97y² = 51 has no solutions Legendre symbol (51/97) = -1 hence 51 is not a quadratic residue of 97.
But this necessary condition is not sufficient...

### Example :

Let equation x² - 7y² = 9
x² - 7y² = 1 has the fundamental solution x0 = 8, y0 = 3 (practice exercice !)
We have to search for 0 ≤ y ≤ 3√(9/18) that is 0 ≤ y ≤ 2.
y = 0 gives x = ±3
y = 1 gives x = ±4
y = 2 gives x² = 37 impossible.

Hence there are 4 primitive solutions, but not all are independant
(±3, 0) are equivalent, more generally (x,y) ~ (-x,-y)
(4,±1) are independant : 16 + 7 = 26 coprime with 9.
(3,0) and (4,1) are independant
(3,0) and (4,-1) are independant

Finally three independant primitive solutions : (3, 0), (4, 1) and (4,-1)
If we just care of x,y > 0, (4,-1) ~ (11,4), is then smallest >0 solution in this family
All solutions are built from the three families :

xn + yn√7 = (8 + 3√7)n(3 + 0√7)
xn + yn√7 = (8 + 3√7)n(4 + √7)
xn + yn√7 = (8 + 3√7)n(4 - √7)

Or the recurrence relations :

xn+1 = 8xn + 21yn,  yn+1 = 3xn + 8yn,    from (x, y) = (3,0), (4,1) and (4,-1)

that is the three families :
(3,0) (24,9) (381,144) ...
(11,4) (172,65) (2741,1036) ...
(4,1) (53,20) (844,319) ...

Scanning for solutions becomes tedious if the bound is too big. In that case we have two efficient methods.

### Case when K² < D

Looking at algorithm for solving x² - Dy² = ±1 again, we have the loop invariant

Pi² - DQi² = (-1)i+1Vi+1V0

This allows to find, with the same algorithm, the primitive solutions of x² - Dy² = K
First of all note that if (x,y) is solution of x² - Dy² = K, (mx,my) is a solution of x² - Dy² = m²K.
If we get a (-1)i+1Vi+1 = K/m², then x = mPi, y = mQi is solution of x² - Dy² = K
If K² < D, this gives all primitive solutions.

Example x² - 19y² = -3

 U V a P Q 0 1 4 4 1 1 4 3 2 9 2 -3 4 1 is solution of x² - 19y² = -3 2 5 1 13 3 5 3 2 3 48 11 -2 3 5 1 61 14 5 2 3 2 170 39 -3 61 14 is solution of x² - 19y² = -3 4 1 8 1 170 39 is solution of x² - 19y² = 1

The primitive solutions of x² - 19y² = -3 are then (4,1) and (61,14)

Notes :
1) Equation x² - 19y² = K hence has no solutions for K = -1, +2, +3, and -4.
2) x² - 19y² = 5 has primitive solutions (9,2) and (48,11) but we don't know if there are others because 5² > 19.

### LMM Algorithm

The Lagrange-Matthews-Mollin algorithm is usefull when K² > D

 List the f  with f² dividing K, m=K/f² For each such f, find all z : -|m|/2 < z≤|m|/2 and z² ≡ D modulo |m| For each z, apply the PQA algorithm with U0 = z, V0 = |m| up to end of period or a Vi = ±1   - If Vi = 1, then (Pi-1,Qi-1) is solution of x² - Dy² = m, add this solution to the list.   - If Vi = -1, then (Pi-1,Qi-1) is solution of x² - Dy² = -m. If x² - Dy² = -1 has a solution (r,s) (primitive)     then (Pi-1,Qi-1)o(r,s) is solution of x² - Dy²=m. Add (Pi-1,Qi-1)o(r,s) = (Pi-1r + Qi-1sD, Pi-1s + Qi-1r)

### Example :

x² - 13y² = 12
f = 1, m = 12, z² ≡ 13 modulo 12, z = ±1 and z = ±5
z = -1 gives the solution (83,23)
z = 1 gives the solution (47,13)
z = -5 gives the solution (5,1)
z = 5 gives the solution (-5,1)
f = 2, m = 3, z² ≡ 13 modulo 3, z = ±1
z = -1 gives the solution (8,2)
z = 1 gives the solution (512,142)

All these solutions being independant, hence result into 6 solution families by combining them with the solution (649,180)n of x² - 13y² = 1

## Equation ax² + bxy + cy² = K

The LMM algorithm can be extended to the more general case of equations ax² + bxy + cy² = K with D = b² - 4ac > 0 not a square.