Calculating the nearest square continued fraction of a quadratic irrational

This program finds the nearest square continued fraction expansion of a quadratic irrational α = (u + t√d)/v, where d,t,u,v are integers, d >1, and nonsquare, t and v nonzero.
We first convert α to (P+√d)/Q, where Q divides d - P2.
Our algorithm is based on papers by A.A. Krisnaswami Ayyangar who defined a reduced surd on page 27 as "the successor of the successor of a special surd".
We use a more direct equivalent version which allows us to detect an NSCF-reduced complete quotient directly.
Let (P + √d)/Q = ξ0 be a standard surd (ie., Q divides d - P2) and gcd(P, Q, (d-P2)/Q)=1). Let a be the integer part of ξ0. Then we can write

ξ0 = a + Q'/(P + √d) or ξ0 = a + 1- Q"/(P" + √d),

where (P' + √d)/Q' and (P" + √d)/Q" are also standard surds.
We choose a0 to be a or a + 1, according as (i) |Q'| < |Q"| or |Q'| > |Q"|, if |Q'| ≠ |Q"| ; or (ii) if |Q'| = |Q"|, according as Q < 0 or Q > 0.

(This differs slightly in case (ii) from Solving the Pell equation, page 407, H.C. Williams, Proc. Millennial Conference on Number Theory, A.K. Peters, 2002, pp. 397-435.)

(We see that |P'2 - d|/|P"2 - d|=|Q'|/|Q"|. Hence |Q'| is ≤ or > |Q"| according as |P'2 - d| is ≤ or > |P"2 - d|.
So choosing the lesser of |Q'| and |Q"| is equivalent to choosing the nearest of the two squares P'2 and P"2 to d.)

Then (P + √d)/Q = ξ0 = a0 + ε1/ ξ1 and where |ε1| = 1 and ξ1 = (P1 + √d)/ Q1 > 1.

We proceed similarly with ξ1 and so on. Then we get complete quotients ξn, where

ξn = an + εn + 1/ ξn + 1 and ξ0 = a0 + ε1/a1 + ε2/a2 + · · ·

We have identities

Pn+1 + Pn = anQn,
Pn+12 + εn+1QnQn+1 = d.

The convergents An/Bn are defined by A-2 = 0, A-1 = 1, B-2 = 1, B-1 = 0 and for k ≥ -1,

Ak+1 = ak+1Ak + εk+1Ak-1,
Bk+1 = ak+1Bk + εk+1Bk-1.

We print the partial numerators (Teilzähler ) εn, the partial denominators (Teilnenner ) an, the convergents and complete quotients, with the period elements in bold font.
Enter d (1 < d < 1010 and not a square):
Enter t (nonzero):
Enter u:
Enter v: (nonzero)

Last modified 29th July 2010
Return to main page