### Solving the Pell equation using NICF

This is a program for finding the least solution of diophantine equations x2 - dy2 = ±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 x2-dy2=-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 q0 = [θ], the nearest integer to θ. i.e. q0 = ⌊θ + 1/2⌋.
Define θk+1 = 1/(qk - θk) and qk+1 = [θk+1].
Then θ = q0 - 1/q1 - 1/q2 - ... - 1/θn and we write θ = (q0, q1, ···, θn) and (q0, q1, ···, qn) = An/Bn.
Then √d has a periodic NICF, with period s + 1:

√d = (q0; q1, ···, qs, q0 + b),
where b = q0 if √d > q0 - γ and b = q0 - 1 , if √d < q0 - γ and γ = (3 - √5)/2 ≈ ·382.

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

Write θk = (Pk + √D)/Qk, so that P0 = 0 and Q0 = 1. Then

qk = [(Pk + √D)/Qk],
Pk+1 = qkQk - Pk,
Qk+1 = (P2k+1 - d)/ Qk.

Also define A-2 = 0, A-1 = 1, B-2 = -1, B-1 = 0 and for k ≥ -1, define

Ak+1 = qk+1Ak - Ak-1
Bk+1 = qk+1Bk - Bk-1.

Then

An/Bn = q0 - 1/q0 - 1/q2 - ... - 1/qn,
Ak-1Bk - Bk-1Ak = 1,
A2k - dB2k = Qk+1.

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

√d = (q0, q1, ···, qm, -q1, -q2, ···, -qm>).

Also (Am-1, Bm-1) is the least solution of x2 - dy2 = -1, while (As-1, Bs-1) is the least solution of x2 - dy2 = 1.

If x2 - dy2 = -1 is not solvable in integers x, y, then s = 2m and again (As-1, Bs-1) is the least solution of x2 - dy2 = 1.

Enter d (< 1010):
Enter e (0 or 1):