Calculating denominators of convergents of a continued fraction

Question

I've written the following Python 2 code to calculate the denominator of the $n$^th continued fraction convergent of $\sqrt{801}$ (this is A042545 on OEIS):

i = input()
a = 801 ** .5
q = [0, 1]
for n in range(i):
  a = 1 / (a % 1)
  q += [int(a) * q[-1] + q[-2]]
print q[i]

This supposedly follows this recursive formula:

$$ q_n = a_n q_{n-1} + q_{n-2}$$

where $a_n$ is calculated as $1$ over the decimal part of $a_{n-1}$. I figured this out myself while studying the sequence, but found confirmation on Wolfram MathWorld (equation 28 on the page).

The program works fine for any input up to $14$, but for input $15$, it produces $27666361$ instead of $53000053$. Interestingly, $16$ produces $53000053$, but the next few results don't seem to match anything in the sequence. I checked WolframAlpha to see if the OEIS entry was possibly wrong, but to no avail.

So either my Python program is incorrect/breaks down at some point, or WolframAlpha is incorrect. The former option seems far more likely, so where am I going wrong?

Calculate first $n$ terms on Ideone. Full output for $30$:

[0, 1, 3, 10, 43, 53, 308, 669, 6998, 7667, 53000, 325667, 2007002, 2332669, 25333692, 27666361, 53000053, 133666467, 1523331190, 3180328847L, 4703660037L, 7883988884L, 12587648921L, 20471637805L, 53530924531L, 181064411398L, 234595335929L, 884850419185L, 1119445755114L, 2004296174299L, 7132334278011L]

Answer 1

Well, here are the correct "digits" $a_n,$ the absolute values of my "deltas" after the initial $28,$ which is not shown by my program. This method, due to Lagrange and Gauss, makes no numerical errors at all, and uses only integer arithmetic. Looking at what you have, it appears you have built some numerical accuracy problems into your calculation.

jagy@phobeusjunior:~/old drive/home/jagy/Cplusplus$ ./Pell 801


0  form   1 56 -17   delta  -3
1  form   -17 46 16   delta  3
2  form   16 50 -11   delta  -4
3  form   -11 38 40   delta  1
4  form   40 42 -9   delta  -5
5  form   -9 48 25   delta  2
6  form   25 52 -5   delta  -10
7  form   -5 48 45   delta  1
8  form   45 42 -8   delta  -6
9  form   -8 54 9   delta  6
10  form   9 54 -8   delta  -6
11  form   -8 42 45   delta  1
12  form   45 48 -5   delta  -10
13  form   -5 52 25   delta  2
14  form   25 48 -9   delta  -5
15  form   -9 42 40   delta  1
16  form   40 38 -11   delta  -4
17  form   -11 50 16   delta  3
18  form   16 46 -17   delta  -3
19  form   -17 56 1   delta  56
20  form   1 56 -17

 disc 3204
Automorph, written on right of Gram matrix:  
5334344001  300333934000
17666702000  994669656001


 Pell automorph 
500002000001  14151028302000
17666702000  500002000001

Pell unit 
500002000001^2 - 801 * 17666702000^2 = 1 

=========================================

801      3^2 *  89

jagy@phobeusjunior:~/old drive/home/jagy/Cplusplus$