### Solving the Pell equation using NICF

This is a program for finding the least solution of diophantine equations
x^{2} - dy^{2} = ±1.

e=1 prints the complete quotients (P+√d)/Q, convergents and partial quotients of the period (or semi-period in the case of even period if x^{2}-dy^{2}=-1 is soluble in integers) of the *nearest integer continued fraction* (NICF) expansion of √d.

Our account is based on exercise 4, page 85 of *Quadratics*, R.A. Mollin, CRC Press 1995

and a nearest-integer formula. Also see *Calculation of the Regulator of Q(√D) by Use of the Nearest Integer Continued Fraction Algorithm*, H.C. Williams and P.A. Buhr, Math. Comp. **33** (1979) 369-381. (We do not use their mid-period improvement.)

Let θ_{0} = θ = √d and q_{0} = [θ], the nearest integer to θ. i.e. q_{0} = ⌊θ + 1/2⌋.

Define &theta_{k+1} = 1/(q_{k} - θ_{k}) and q_{k+1} = [&theta_{k+1}].

Then θ = q_{0} - 1/q_{1} - 1/q_{2} - ... - 1/θ_{n} and we write θ
= (q_{0}, q_{1}, ···, θ_{n}) and (q_{0}, q_{1}, ···, q_{n}) = A_{n}/B_{n}.

Then √d has a periodic NICF, with period s + 1:

√d = (q_{0}; q_{1}, ···, q_{s}, q_{0} + b),

where b = q_{0} if √d > q_{0} - γ and b = q_{0} - 1 , if √d < q_{0} - γ and γ = (3 - √5)/2 ≈ ·382.

d = 13 is the first example where b = q_{0} - 1: √13 = (4; 3, 2, -7, -3, -2, 7).
Write &theta_{k} = (P_{k} + √D)/Q_{k}, so that
P_{0} = 0 and Q_{0} = 1. Then

q_{k} = [(P_{k} + √D)/Q_{k}],

P_{k+1} = q_{k}Q_{k} - P_{k},

Q_{k+1} = (P^{2}_{k+1} - d)/ Q_{k}.

Also define
A_{-2} = 0,
A_{-1} = 1,
B_{-2} = -1,
B_{-1} = 0 and for k ≥ -1, define
A_{k+1} = q_{k+1}A_{k} - A_{k-1}

B_{k+1} = q_{k+1}B_{k} - B_{k-1}.

Then
A_{n}/B_{n} = q_{0} - 1/q_{0} - 1/q_{2} - ... - 1/q_{n},

A_{k-1}B_{k} - B_{k-1}A_{k} = 1,

A^{2}_{k} - dB^{2}_{k} = Q_{k+1}.

If x^{2} - dy^{2} = -1 is solvable in integers x, y, then s = 2m - 1 and the NICF of √d is
√d = (q_{0}, \overline{q_{1}, ···, q_{m}, -q_{1}, -q_{2}, ···, -q_{m}}).

Also (A_{m-1}, B_{m-1}) is the least solution of x^{2} - dy^{2} = -1, while (A_{s-1}, B_{s-1}) is the least solution of x^{2} - dy^{2} = 1.
If x^{2} - dy^{2} = -1 is not solvable in integers x, y, then s = 2m and again (A_{s-1}, B_{s-1}) is the least solution of x^{2} - dy^{2} = 1.

*Last modified 7th August 2008*

Return to main page