﻿ Euler Math Toolbox - 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 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
```

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

Examples