#### Calculating a continued fraction arising from an exceptional solution to Andrej Dujella's diophantine equation x2 − (k2+1)y2 = k2

We calculate the continued fraction of √(b/a), where b > a > 2, ab = k2 + 1, gcd(a, b) = 1.
Here a and b correspond to an exceptional solution (x, y) of Dujella's equation, i.e., where 1 ≤ y < k − 1:
if d = gcd(x + k, x − k), then a = gcd((x + k)/d, k2 + 1), b = gcd((x − k)/d, k2 + 1). See slides.

If a > b, it is convenient to interchange a and b.

All exceptional solutions are given by a recursively defined forest, starting from three types of root nodes:

• (t,t,0), t ≥ 2 (type 0),
• (t, t2 + t + 1, t + 1), t ≥ 1 (type 1)
• (t, t2 − t + 1, t − 1), t ≥ 2 (type -1).
A ternary expansion of n (digits 1, 0 or -1) determines a branch of the forest giving an exceptional solution as follows:
We find n = r 0 + 3r1 + ··· + 3u −1ru −1 + 3u.
Then the exceptional solution (k,x,y) is obtained from a root node by applying first g+, then successively g+, g0 or g-, according as ri = 1, 0 or -1, for i = u-1,...,0.
We print the digit string r0r1···ru-11.

It is classical that the continued fraction is periodic after the first partial quotient. Here the period-length = 2h + 1, with partial quotients a0, a1,..., ah, ah,..., a1, 2a0,... We print only a0, a1,..., ah.

This is done for t, t+1, t+2 and t+3, so as to allow guessing the linear polynomials in t that conjecturally form the partial quotients. The polynomial partial quotients are also given.

Three phenomena are clear experimentally:

1. For exceptional solutions (k(t),x(t),y(t)) with root node (t,t,0), the period length is constant for t ≥ 2.
2. For exceptional solutions (k(t),x(t),y(t)) with root node (t,t2+et+1,t+e), e=±1, the period length for d = 1 is constant, say L, for t ≥ 2 and parity of t, while that for d = 2 is equal to 3L for t > 2 and opposite parity.
3. If d = 1, all partial quotients are even.
Enter t (≥ 2 if type = 0 or -1, ≥ 1 if type = 1):
Enter n (1 ≤ 10000):
Enter type (1,0 or -1 ):