Computing continued fraction expansions

Question

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).

Answer 1

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 $$

Answer 2

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)

Answer 3

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
=========================================

Answer 4

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]