# Computing continued fraction expansions

My question concerns the numerical accuracy of a continued fraction expansion. A typical algorithm for computing a continued fraction can be written in Python as :

x0 = sqrt(2)

N = 40
a = [0]*N
u = [0]*N

x = x0
for k in range(N):
a[k] = int(x//1)
u[k] = x % 1      # Often replaced with x - a[k]  ???
x = 1/u[k]

print(a)


For $$x=\sqrt{2}$$, this produces

[1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 3, 3, 1, 3, 1, 1, 2, 1809, 1, 2, 5, 2, 2, 1, 2, 1]


Or for $$x = \pi$$, I get :

[3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 3, 3, 23, 1, 1, 7, 4, 35, 1, 1, 1, 2, 3, 3, 3, 3, 1, 1, 14, 6, 4, 5, 1, 7, 1, 5, 1]


One observation is that round-off error has clearly crept into the calculation. The expansion for $$\sqrt{2}$$ should be repeating 2s, and the approximation for $$\pi$$ is accurate up to the first appearance of 14 (12 terms after the leading $$a_0=3$$ term).

What I also found is that while round off error does creep in, the number of correct terms in the expansion is exactly what is needed to get convergents that agree with $$x$$ up to machine precision. So for example, convergent 20 of $$\sqrt{2}$$ is given by

54608393 / 38613965


This convergent approximates $$\sqrt{2}$$ to $$4.44 \times 10^{-16}$$. Convergent 21 agrees with $$\sqrt{2}$$ exactly (in finite precision arithmetic). Interestingly, the expansion has exactly 20 correct terms (i.e. 2s) after the leading $$a_0 = 1$$.

Similarly, for $$\pi$$, the convergent 12 is given by

80143857 / 25510582


which agrees with $$\pi$$ to within an error of $$4.44 \times 10^{-16}$$. Convergent 13 agrees with $$\pi$$ exactly (in finite precision arithmetic).

In both cases, the numerators for the convergents agree with the OEIS. See A001333 and A002485.

The above observations led me to these questions :

• Is there a standard stopping criteria that can be used to determine when the expansion has converged to a desired tolerance?

• Is it always the case that one will have enough correct terms in the expansion to approximate the desired number to machine precision?

• Is it possible to detect whether a continued fraction is periodic? Or, if one knows it will be periodic (i.e. $$x$$ is a root of a quadratic with integer coefficients), is it possible to get the periodic sequence exactly?

It has also occurred to me that nobody would think of computing continued fraction expansions using finite precision arithmetic! I would like to do this problem as an exercise in a beginning computational math course, and was hoping not to go into variable precision arithmetic (an area which is really outside of my wheel house).

• there are methods for the (simple) continued fraction for $\sqrt n$ and for $\frac{a+\sqrt b}{c}$ that require only integer operations and produce no errors. For anything more complicated, it is a matter of luck whether there is any nice way to get perfect accuracy; Commented Mar 1, 2021 at 2:49
• Let's see, you do ask about periodic; a purely periodic fraction converges to a quadratic irrational that is also reduced. For example, $\lfloor \sqrt n \rfloor + \sqrt n$ for positive integer $n$ is reduced (here $n$ is not a square) Commented Mar 1, 2021 at 2:54
• the most general statements about error are in Khinchin's little book. here is one by Olds ms.uky.edu/~sohum/ma330/files/Continued%20Fractions.pdf Commented Mar 1, 2021 at 17:02
• math.ru.nl/~bosma/Students/CF.pdf Commented Mar 1, 2021 at 17:04
• amazon.com/Continued-Fractions-Dover-Books-Mathematics/dp/… Commented Mar 1, 2021 at 17:07

Here is a method for $$\sqrt n$$ that is likely to be what Fermat used in hand computations.

$$\sqrt { 5} = 2 + \frac{ \sqrt {5} - 2 }{ 1 }$$ $$\frac{ 1 }{ \sqrt {5} - 2 } = \frac{ \sqrt {5} + 2 }{1 } = 4 + \frac{ \sqrt {5} - 2 }{1 }$$

Simple continued fraction tableau:
$$\begin{array}{cccccccc} & & 2 & & 4 & & 4 & \\ \\ \frac{ 0 }{ 1 } & \frac{ 1 }{ 0 } & & \frac{ 2 }{ 1 } & & \frac{ 9 }{ 4 } \\ \\ & 1 & & -1 & & 1 \end{array}$$

$$\begin{array}{cccc} \frac{ 1 }{ 0 } & 1^2 - 5 \cdot 0^2 = 1 & \mbox{digit} & 2 \\ \frac{ 2 }{ 1 } & 2^2 - 5 \cdot 1^2 = -1 & \mbox{digit} & 4 \\ \frac{ 9 }{ 4 } & 9^2 - 5 \cdot 4^2 = 1 & \mbox{digit} & 4 \\ \end{array}$$

$$\bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc\bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc$$

$$\sqrt { 13} = 3 + \frac{ \sqrt {13} - 3 }{ 1 }$$ $$\frac{ 1 }{ \sqrt {13} - 3 } = \frac{ \sqrt {13} + 3 }{4 } = 1 + \frac{ \sqrt {13} - 1 }{4 }$$ $$\frac{ 4 }{ \sqrt {13} - 1 } = \frac{ \sqrt {13} + 1 }{3 } = 1 + \frac{ \sqrt {13} - 2 }{3 }$$ $$\frac{ 3 }{ \sqrt {13} - 2 } = \frac{ \sqrt {13} + 2 }{3 } = 1 + \frac{ \sqrt {13} - 1 }{3 }$$ $$\frac{ 3 }{ \sqrt {13} - 1 } = \frac{ \sqrt {13} + 1 }{4 } = 1 + \frac{ \sqrt {13} - 3 }{4 }$$ $$\frac{ 4 }{ \sqrt {13} - 3 } = \frac{ \sqrt {13} + 3 }{1 } = 6 + \frac{ \sqrt {13} - 3 }{1 }$$

Simple continued fraction tableau:
$$\begin{array}{cccccccccccccccccccccccc} & & 3 & & 1 & & 1 & & 1 & & 1 & & 6 & & 1 & & 1 & & 1 & & 1 & & 6 & \\ \\ \frac{ 0 }{ 1 } & \frac{ 1 }{ 0 } & & \frac{ 3 }{ 1 } & & \frac{ 4 }{ 1 } & & \frac{ 7 }{ 2 } & & \frac{ 11 }{ 3 } & & \frac{ 18 }{ 5 } & & \frac{ 119 }{ 33 } & & \frac{ 137 }{ 38 } & & \frac{ 256 }{ 71 } & & \frac{ 393 }{ 109 } & & \frac{ 649 }{ 180 } \\ \\ & 1 & & -4 & & 3 & & -3 & & 4 & & -1 & & 4 & & -3 & & 3 & & -4 & & 1 \end{array}$$

$$\begin{array}{cccc} \frac{ 1 }{ 0 } & 1^2 - 13 \cdot 0^2 = 1 & \mbox{digit} & 3 \\ \frac{ 3 }{ 1 } & 3^2 - 13 \cdot 1^2 = -4 & \mbox{digit} & 1 \\ \frac{ 4 }{ 1 } & 4^2 - 13 \cdot 1^2 = 3 & \mbox{digit} & 1 \\ \frac{ 7 }{ 2 } & 7^2 - 13 \cdot 2^2 = -3 & \mbox{digit} & 1 \\ \frac{ 11 }{ 3 } & 11^2 - 13 \cdot 3^2 = 4 & \mbox{digit} & 1 \\ \frac{ 18 }{ 5 } & 18^2 - 13 \cdot 5^2 = -1 & \mbox{digit} & 6 \\ \frac{ 119 }{ 33 } & 119^2 - 13 \cdot 33^2 = 4 & \mbox{digit} & 1 \\ \frac{ 137 }{ 38 } & 137^2 - 13 \cdot 38^2 = -3 & \mbox{digit} & 1 \\ \frac{ 256 }{ 71 } & 256^2 - 13 \cdot 71^2 = 3 & \mbox{digit} & 1 \\ \frac{ 393 }{ 109 } & 393^2 - 13 \cdot 109^2 = -4 & \mbox{digit} & 1 \\ \frac{ 649 }{ 180 } & 649^2 - 13 \cdot 180^2 = 1 & \mbox{digit} & 6 \\ \end{array}$$

$$\bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc\bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc$$

• I get the first two steps in each equality (just rationalizing the denominator) But how are you choosing to break up the single fraction into an integer plus a fraction? I see that the goal is to get back to a term in the original "add-subtract" step (sqrt(13)-3)/1 in your second example. Commented Mar 6, 2021 at 17:10
• Method described by Prof. Lubin at math.stackexchange.com/questions/2215918/… Commented Mar 6, 2021 at 17:18

Here is a program and the protocol, using Pari/GP; but I hope the program is understandable. The stopping-criterion is, whether the error in the reconstruction of the intended irrational number is equal the machine-epsilon.
I set the Pari/GP precision to small number of digits and then run the loop. Using $$\sqrt2$$.

fmt(10)       \\ set very small ;-) internal precision (19 dec digits)
va=vector(50) \\ this takes the found coefficients, set length as you think serves well
w=b=sqrt(2)   \\ this is the initial value, shall have 18 or 19 dec digits correct
ia=0          \\ through iterating, this is the index for the va-vector
pq=[0,1;1,0]  \\ this captures the current convergents (p,q)
\\  allows checking whether  w - p/q  <= machine-epsilon



protocol:

\\  w=  1.41421356237 \\ protocol of w=b=sqrt(2), internal 18 or 19 dec digits correct


program:

{for(k=1,50,  \\ set upper limit as you think serves well
a=floor(b);
ia++; va[ia]=a; \\ save the cf-digit (="cf-quotient") in vector
b=1/(b-a);
pq=pq*[0,1;1,a];  \\ compute current convergent
p=pq[1,2];q=pq[2,2];
print([k,a,err = w - p/q]); \\ document progress
if(err==0.0,break());   \\ if machine-epsilon reached, stop iteration
);}



Protocol:

 k  a  err             :: a=current quotient/digit, err=w - p/q
[1, 1, 0.414213562373]
[2, 2, -0.0857864376269]
[3, 2, 0.0142135623731]
[4, 2, -0.00245310429357]
[5, 2, 0.000420458924819]
[6, 2, -0.0000721519126192]
[7, 2, 0.0000123789411424]
[8, 2, -0.00000212390141476]
[9, 2, 0.000000364403551901]
[10, 2, -0.0000000625217745896]
[11, 2, 0.0000000107270403544]
[12, 2, -0.00000000184046916492]
[13, 2, 0.0000000003157745862]
[14, 2, -5.417835374 E-11]
[15, 2, 9.29553586 E-12]
[16, 2, -1.594861909 E-12]
[17, 2, 2.736349747 E-13]
[18, 2, -4.694840810 E-14]
[19, 2, 8.05490168 E-15]
[20, 2, -1.382117461 E-15]
[21, 2, 2.370426006 E-16]
[22, 2, -4.074183443 E-17]
[23, 2, 6.853781456 E-18]
[24, 2, -1.285552860 E-18]
[25, 2, 8.87489049 E-20]
[26, 2, -1.470440003 E-19]
[27, 1, 0.E-18]             \\ stopped at iteration 27


Show continued-fraction quotients

 va=vecextract(va,Str("1..",ia))} \\ remove superfluous trailing entries in va
%630 = [1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 1]
-----------------------------------^ last entry is likely wrong


Using $$w=\pi$$ gives the following protocol

[1, 3, 0.141592653590]
[2, 7, -0.00126448926735]
[3, 15, 0.0000832196275291]
[4, 1, -0.000000266764189063]
[5, 292, 0.000000000577890634]
[6, 1, -0.0000000003316278063]
[7, 1, 1.223565329 E-10]
[8, 1, -2.914338504 E-11]
[9, 2, 8.71546711 E-12]
[10, 1, -1.610740232 E-12]
[11, 3, 4.040668867 E-13]
[12, 1, -2.214484967 E-14]
[13, 14, 5.78983025 E-16]
[14, 2, -1.641519209 E-16]
[15, 1, 7.80993605 E-17]
[16, 1, -1.932805900 E-17]
[17, 2, 3.079202631 E-18]
[18, 2, -5.73807725 E-19]
[19, 2, 0.E-18]


and the continued fraction

 %645 = [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2]
------------------------------------------------------------------^
(last quotient a=2 likely wrong)

• This looks like the best approach. Initially, I thought that it was possible to terminate the sequence just by knowing the sequence a. But realize now that using the convergents is the obvious way to go. Your computation of the convergents via matrices (I think) is nice. It would still be an interesting exercise to understand where floating point arithmetic fails in the computation of the expansion. I guess it is in the repeated truncation and inverting, until only garbage digits are left. But 16 will always be good. Commented Mar 5, 2021 at 22:46
• @Donna - Pari/GP seems to be reliable in this matter. Just tried my old languages Turbo-Pascal and Delphi. In Turbo-Pascal the implementation of float values allows the "epsilon = 0.0" - test, but only if the float-datatype is "real". If not, the eps=0.0-case need not occur at all! (and garbage-cf-coefficients will be produced) The test for last valid quotient/digit must then include, whether the current absolute value of error is not larger than previous absolute error. It seems this suffices for the criterion to leave the iteration-loop with last quotient/digit insecure. Commented Mar 6, 2021 at 0:11
• It is interesting that "err == 0" works at all. This would normally be a rookie mistake in numerical analysis (we never get exactly zero!) But here, it seems that not only do the convergents hit the 64-bit representation exactly, they do not get any worse if one continues the iteration beyond this point. Commented Mar 6, 2021 at 15:58
• @Donna - just a nitpick. In Pari/GP "0" and "0.0" is not equal. If I had written if(err==0,...) this would likely not work. I think, that Pari/GP actively "watches" the 19-digit internal precision (which I demanded in the beginning) over the successive operations and "sees" the exhausting of significant digits- I've no other explanation so far. But one could ask the Pari/GP-community for more exact/internal information... Commented Mar 6, 2021 at 20:51
• Thanks - I didn't pick up on that (not being entirely sure what Pari/GP is). But at least in Python, the error in the convergents really does seem to hit exactly 0. Maybe this isn't too surprising, though. Commented Mar 7, 2021 at 22:48

This is the Gauss-Lagrange method of neighboring reduced indefinite forms. I learned this in Binary Quadratic Forms by D. A. Buell. It is also in Dickson (1929) Introduction to the Theory of Numbers.

The output below says that a root of $$7x^2 + 3 xy - 7 y^2$$ has purely periodic fraction 1,4,4,1. (You just take the absolute values of my "delta" numbers to make the continued fraction). Let me think about which root that might be. Alright, got it. It is actually the ratio $$r=y/x$$ which solves $$7+3r-7r^2,$$ the relevant root being $$\frac{3+ \sqrt{205}}{14}$$

=========================================
jagy@phobeusjunior:~/old drive/home/jagy/Cplusplus$./indefCycle 7 3 -7 0 form 7 3 -7 1 0 0 1 To Return 1 0 0 1 0 form 7 3 -7 delta -1 1 form -7 11 3 delta 4 2 form 3 13 -3 delta -4 3 form -3 11 7 delta 1 4 form 7 3 -7 form 7 x^2 + 3 x y -7 y^2 minimum was 3rep x = -1 y = -1 disc 205 dSqrt 14 M_Ratio 4 Automorph, written on right of Gram matrix: 17 21 21 26 =========================================  This one is for $$\frac{13 + \sqrt{1313}}{44}$$ jagy@phobeusjunior:~/old drive/home/jagy/Cplusplus$ ./indefCycle 13 13 -22

0  form             13          13         -22

1           0
0           1

To Return
1           0
0           1

0  form   13 13 -22   delta  -1     ambiguous
1  form   -22 31 4   delta  8
2  form   4 33 -14   delta  -2
3  form   -14 23 14   delta  2
4  form   14 33 -4   delta  -8
5  form   -4 31 22   delta  1
6  form   22 13 -13   delta  -1
7  form   -13 13 22   delta  1     ambiguous            -1 composed with form zero
8  form   22 31 -4   delta  -8
9  form   -4 33 14   delta  2
10  form   14 23 -14   delta  -2
11  form   -14 33 4   delta  8
12  form   4 31 -22   delta  -1
13  form   -22 13 13   delta  1
14  form   13 13 -22

form   13 x^2  + 13 x y  -22 y^2

minimum was   4rep   x = -1   y = -1 disc 1313 dSqrt 36  M_Ratio  7.668639
Automorph, written on right of Gram matrix:
-486641  -921536
-544544  -1031185
=========================================


took a while, I wrote the Gauss-Lagrange method using the left neighbor instead of the right neighbor. Given a reduced form $$Ax^2 + B xy + C y^2,$$ which means that $$D= B^2 - 4 AC$$ is positive but not a square, $$AC < 0$$ and $$B > |A+C|,$$ the numbers listed that I call "epsilon" (well, its absolute value) give the partial quotients for the continued fraction of $$\frac{B+ \sqrt D}{2A}$$ when $$A$$ is positive. I have not tried negative $$A$$ yet.. sample run

Sample: $$47x^2 + 321 xy - 37 y^2$$ and $$D=109997,$$ target is $$\frac{321+ \sqrt{ 109997}}{94}$$

First, a few "partial quotients" in the continued fraction expansion

? x = 321 + sqrt(109997)
%3 = 652.6579563345345104401576022
? x /= 94
%4 = 6.943169748239728834469761725
?
? x -= floor(x) ; x = 1/x
%5 = 1.060254531982535997122061997
? x -= floor(x) ; x = 1/x
%6 = 16.59626200880353974842520003
? x -= floor(x) ; x = 1/x
%7 = 1.677115068938572016859051397
? x -= floor(x) ; x = 1/x
%8 = 1.476853855235527410365538515
? x -= floor(x) ; x = 1/x
%9 = 2.097078568246198860000621270
? x -= floor(x) ; x = 1/x
%10 = 10.30093477958926629742181316
? x -= floor(x) ; x = 1/x
%11 = 3.322979156363579950723385592

jagy@phobeusjunior:~/old drive/home/jagy/Cplusplus\$ ./indefCycleLeft  47  321  -37

0  form   47 321 -37   epsilon  6
1  form   -271 243 47   Epsilon  -1
2  form   19 299 -271   Epsilon  16
3  form   -191 309 19   Epsilon  -1
4  form   137 73 -191   Epsilon  1
5  form   -127 201 137   Epsilon  -2
6  form   31 307 -127   Epsilon  10
7  form   -97 313 31   Epsilon  -3
8  form   97 269 -97   Epsilon  3
9  form   -31 313 97   Epsilon  -10
10  form   127 307 -31   Epsilon  2
11  form   -137 201 127   Epsilon  -1
12  form   191 73 -137   Epsilon  1
13  form   -19 309 191   Epsilon  -16
14  form   271 299 -19   Epsilon  1
15  form   -47 243 271   Epsilon  -6
16  form   37 321 -47   Epsilon  8
17  form   -247 271 37   Epsilon  -1
18  form   61 223 -247   Epsilon  4
19  form   -163 265 61   Epsilon  -1
20  form   163 61 -163   Epsilon  1
21  form   -61 265 163   Epsilon  -4
22  form   247 223 -61   Epsilon  1
23  form   -37 271 247   Epsilon  -8
24  form   47 321 -37
form   47 x^2  + 321 x y  -37 y^2

minimum was   19rep   x = -7   y = 1 disc 109997 dSqrt 331
Automorph, written on right of Gram matrix:
2060501320695  -233624820247
-296766663557  33648150444
for  gp Pari: rt =  [ 2060501320695 , -296766663557 ; -233624820247 , 33648150444 ] ;    h =  [ 94 , 321 ; 321 , -74 ] ;    r =  [ 2060501320695 , -233624820247 ; -296766663557 , 33648150444 ] ;
=========================================

PARI/GP is free software, covered by the GNU General Public License,
?
? rt
%7 =
[2060501320695 -296766663557]
[-233624820247   33648150444]

? h
%8 =
[ 94 321]
[321 -74]

? r
%9 =
[2060501320695 -233624820247]
[-296766663557   33648150444]

? rt * h * r
%10 =
[ 94 321]
[321 -74]