Solving the Pell equation using NICF

This is a program for finding the least solution of diophantine equations x² - dy² = ±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²-dy²=-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 θ₀ = θ = √d and q₀ = [θ], the nearest integer to θ. i.e. q₀ = ⌊θ + 1/2⌋.
Define θ_k+1 = 1/(q_k - θ_k) and q_k+1 = [θ_k+1].
Then θ = q₀ - 1/q₁ - 1/q₂ - ... - 1/θ_n and we write θ = (q₀, q₁, ···, θ_n) and (q₀, q₁, ···, q_n) = A_n/B_n.
Then √d has a periodic NICF, with period s + 1:

√d = (q₀; q₁, ···, q_s, q₀ + b),
where b = q₀ if √d > q₀ - γ and b = q₀ - 1 , if √d < q₀ - γ and γ = (3 - √5)/2 ≈ ·382.

d = 13 is the first example where b = q₀ - 1: √13 = (4; 3, 2, -7, -3, -2, 7).

Write θ_k = (P_k + √D)/Q_k, so that P₀ = 0 and Q₀ = 1. Then

q_k = [(P_k + √D)/Q_k],
P_k+1 = q_kQ_k - P_k,
Q_k+1 = (P²_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+1A_k - A_k-1
B_k+1 = q_k+1B_k - B_k-1.

Then

A_n/B_n = q₀ - 1/q₀ - 1/q₂ - ... - 1/q_n,
A_k-1B_k - B_k-1A_k = 1,
A²_k - dB²_k = Q_k+1.

If x² - dy² = -1 is solvable in integers x, y, then s = 2m - 1 and the NICF of √d is

√d = (q₀, q₁, ···, q_m, -q₁, -q₂, ···, -q_m>).

Also (A_m-1, B_m-1) is the least solution of x² - dy² = -1, while (A_s-1, B_s-1) is the least solution of x² - dy² = 1.

If x² - dy² = -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² - dy² = 1.

Last modified 7th August 2008
Return to main page