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 function of Euler.

>contfrac(sqrt(2),2)

[1, 2, 2]

The vector represents the fraction

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

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

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. The function fraction prints a real number as a fraction (using continued fractions internally).

>fraction contfracbest(pi,3)

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.

>fraction contfracbest(pi,1)

22/7

>22/7-pi

0.00126448926735

The continued fraction of E has interesting coefficients. We do not prove this here, but it follows from the series of E. It can be used to show that E is transcendental.

>contfrac(E,12)

[2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1]

>fraction contfracbest(E,8)

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[-1],2*n-1)

[1, 1.5, 1.41667] [1, 2, 2, 2]

>n=3; h=iterate("(x+2/x)/2",1,n); contfrac(h[-1],2*n-1)

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

But the length of the approximation grows 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