iconEuler Examples

Continued Fractions

by R. Grothmann

In this notebook, we study continued fractions. A continued fraction
(CF) is a way to approximate a real number with a fraction. E.g.,
ancient cultures used pi=22/7 as a very close approximation. How can
we find such approximations?

The following is the well known CF expansion of sqrt(2) computed with
the contfrac funciton of Euler.
>contfrac(sqrt(2),2)
[ 1  2  2 ]
The vector represents the fraction

Continued Fractions

We can evaluate this, and get two approximating fractions, one with
(+1) and one without.
>frac(1+1/(2+1/2)), frac(1+1/(2+1/3))
7/5
10/7
One of these approximation is the best possible rational approximation
with the same or a smaller numerator.

In this case the first is the better approximation.
>1+1/(2+1/2)-sqrt(2), 1+1/(2+1/3)-sqrt(2)
-0.0142135623731
0.0143578661983
It is not difficult to get the best approximation with a fixed
numerator, since

Continued Fractions

So we only have to round mx to the closest integer. The following
function does this with an effort of O(m) for m<=M.
>function bestrat (x,M) ...
nbest=round(x,0);
mbest=1;
err=abs(x-nbest/mbest);
for m=2 to M;
  n=round(m*x,0);
  if abs(x-n/m)<err then
    nbest=n; mbest=m; err=abs(x-n/m);
  endif;
end;
return {nbest,mbest}
endfunction
Let us test with sqrt(2) and m=7.
>{n,m}=bestrat(sqrt(2),7); n+"/"+m
7/5
The following approximation of pi has been known since ancient times.
>{n,m}=bestrat(pi,500); n+"/"+m
355/113
The algorithm to compute continued fractions works differently. It is
much faster, and it can be generalized for rational functions or other
fields.

The idea is to set

Continued Fractions

Continued Fractions

Then we continue recursively with r_1. floor(x) is the integer part of
x, of course.

The following simple recursion does this, printing the coefficients
a_i along the way.
>function docontfrac (x,n) ...
a=floor(x),
if n>0 then docontfrac(1/(x-a),n-1); endif;
endfunction
>docontfrac(sqrt(2),5)
1
2
2
2
2
2
The built-in contfrac function uses a loop instead. But otherwise it
works in the same way.
>contfrac(sqrt(2),5)
[ 1  2  2  2  2  2 ]
The built-in function contfracval evaluates both approximations. It
uses a loop too. The implementation is simple.
>type contfracval
function contfracval (r)
    n=cols(r);
    x1=r[n]; x2=r[n]+1;
    loop 1 to n-1
         x1=1/x1+r[n-#];
         x2=1/x2+r[n-#];
    end;
    return {x1,x2};
endfunction
It returns two real numbers. We use the frac(x) command to print them
in fractional format (which in turn uses continued fractions).
>{a,b}=contfracval(contfrac(sqrt(2),5)); frac(a), frac(b),
99/70
140/99
The same for pi.
>{a,b}=contfracval(contfrac(pi,3)); frac(a), frac(b),
355/113
688/219
We can automatically choose the better one.
>fracformat; contfracbest(pi,3), longformat;
355/113
The error is good enough for most earthly purposes.
>355/113-pi
2.66764189405e-007
Here is a shorter CF, which is also an ancient approximation.
>fracformat; contfracbest(pi,1), longformat;
22/7
>22/7-pi
0.00126448926735
The continue fraction of E has increasing coefficients.
>contfrac(E,12)
[ 2  1  2  1  1  4  1  1  6  1  1  8  1 ]
>fracformat; contfracbest(E,8), longformat;
1457/536
>1457/536-E
1.75363050703e-006
The CF expansion of sqrt(2) has the following explanation. The well
known Newton iteration to sqrt(2) produces exactly the CF
approximations, when started in 1.
>n=2; h=iterate("(x+2/x)/2",1,n); contfrac(h[n],2*n-1)
[ 1  2  2  2 ]
>n=3; h=iterate("(x+2/x)/2",1,n); contfrac(h[n],2*n-1)
[ 1  2  2  2  2  2 ]
But the length of the approximation grows like exponentially.
>n=5; h=iterate("(x+2/x)/2",1,n); contfrac(h[n],2*n-1)
[ 1  2  2  2  2  2  2  2  2  2 ]
Let us try to compute this process in Maxima.
>function f(x) &= (x+2/x)/2
                                    2
                                x + -
                                    x
                                -----
                                  2

Maxima has an algorithm for continued fractions and its evaluation
too.
>a &= f(f(f(1))), &cf(a), &ratsimp(cfdisrep(%))
                                 577
                                 ---
                                 408


                       [1, 2, 2, 2, 2, 2, 2, 2]


                                 577
                                 ---
                                 408

There is a famous approximation of the halftone interval in music.
>fracprint(contfracbest(2^(1/12),2));
18/17
It is a bit two short, when applied 12 times.
>(18/17)^12
1.98555995207
But the half tone is only 1 cent short, which is very good.
>1200*(log(18/17)-log(2^(1/12)))/log(2)
-1.04540776963
12 such half tones are 10 cent short.
>1200*(12*log(18/17)-log(2))/log(2)
-12.5448932356
This is the next best approximation.
>fracprint(contfracbest(2^(1/12),3));
107/101

Examples Homepage