# Continued Fractions and Function “rat” 1

Posted by **Cleve Moler**,

Let me tell you about MATLAB's controversial function `rat`.

### Contents

#### Classic `rat`

The function `rat` is older than MathWorks. It was one of the 73 functions in my original Fortran MATLAB. Here is the help entry from 40 years ago.

```
<>
help rat
```

RAT An experimental function which attempts to remove the roundoff error from results that should be "simple" rational numbers. RAT(X) approximates each element of X by a continued fraction of the form

a/b = d1 + 1/(d2 + 1/(d3 + ... + 1/dk))

with k <= len, integer di and abs(di) <= max . The default values of the parameters are len = 5 and max = 100. RAT(len,max) changes the default values. Increasing either len or max increases the number of possible fractions. <A,B> = RAT(X) produces integer matrices A and B so that

A ./ B = RAT(X)

Some examples:

long T = hilb(6), X = inv(T) <A,B> = rat(X) H = A ./ B, S = inv(H)

short e d = 1:8, e = ones(d), A = abs(d'*e - e'*d) X = inv(A) rat(X) display(ans)

#### An example

Let's try that second example with modern MATLAB and `format rat`.

format short d = 1:6 e = ones(size(d)), A = abs(d'*e - e'*d) format short e X = inv(A) format rat X

d = 1 2 3 4 5 6 e = 1 1 1 1 1 1 A = 0 1 2 3 4 5 1 0 1 2 3 4 2 1 0 1 2 3 3 2 1 0 1 2 4 3 2 1 0 1 5 4 3 2 1 0 X = -4.0000e-01 5.0000e-01 8.3267e-17 -3.3307e-17 -9.9920e-17 1.0000e-01 5.0000e-01 -1.0000e+00 5.0000e-01 3.6978e-33 8.8818e-17 -5.5511e-17 8.3267e-17 5.0000e-01 -1.0000e+00 5.0000e-01 1.4063e-16 -9.2519e-17 -3.3307e-17 3.6978e-33 5.0000e-01 -1.0000e+00 5.0000e-01 5.5511e-17 -9.9920e-17 8.8818e-17 1.4063e-16 5.0000e-01 -1.0000e+00 5.0000e-01 1.0000e-01 -5.5511e-17 -9.2519e-17 5.5511e-17 5.0000e-01 -4.0000e-01 X = Columns 1 through 5 -2/5 1/2 * * * 1/2 -1 1/2 * * * 1/2 -1 1/2 * * * 1/2 -1 1/2 * * * 1/2 -1 1/10 * * * 1/2 Column 6 1/10 * * * 1/2 -2/5

Which version of `X` do you prefer? The first, printed with a floating point format so the roundoff errors are readily apparent. Or the second, printed with a rational format that seeks to disguise the roundoff. For this example, almost everybody would vote for the rational format.

#### Don't try to hide it

But I have learned that it is better to explain floating point arithmetic than it is to try to hide it. Besides, most important MATLAB calculations involve larger matrices, and eigenvalues, differential equations and signal processing. They don't have such tidy results.

So, we don't have much use for `rat` and `format rat` today, but they're still there and, as I said, they're moderately controversial.

#### The controversial algorithm

The core of `rat` generates continued fractions by repeatedly subtracting off the integer part and taking the reciprocal of what is left. If there is nothing left to reciprocate, the process terminates because the input is a rational number. If the input is irrational, the process does not terminate.

The controversary comes from what is meant by integer part. To generate proper continued fractions, integer part should be `floor`. That always leaves a positive fraction to reciprocate. But 40 years ago, I got clever and used `round` instead of `floor`. That means it may take fewer terms to obtain a specified accuracy, but the continued fractions are, shall I say, *unorthodox*. Let's see more examples.

#### pi

Let's generate increasingly accurate rational approximations to $\pi$. With two output arguments, `rat` unwinds the continued fraction to produce two integers whose ratio has the same value. The accuracy of `rat` is determined by an optional tolerance, `rat(x,tol)`. The default tolerance is `1.0e-6*norm(x)`.

Requesting low accuracy produces a familiar approximation of $\pi$.

```
form = ' tol = %6.1e \n r = %s \n = %d/%d\n';
tol = 1.0e-2;
r = rat(pi,tol);
[n,d] = rat(pi,tol);
fprintf(form,tol,r,n,d)
```

tol = 1.0e-02 r = 3 + 1/(7) = 22/7

The default tolerance produces a less familiar, but elegant approximation.

tol = 1.0e-6; r = rat(pi,tol); [n,d] = rat(pi,tol); fprintf(form,tol,r,n,d)

tol = 1.0e-06 r = 3 + 1/(7 + 1/(16)) = 355/113

Somewhat higher accuracy produces a negative term.

tol = 1.0e-8; r = rat(pi,tol); [n,d] = rat(pi,tol); fprintf(form,tol,r,n,d)

tol = 1.0e-08 r = 3 + 1/(7 + 1/(16 + 1/(-294))) = 104348/33215

Another negative term.

tol = 1.0e-11; r = rat(pi,tol); [n,d] = rat(pi,tol); fprintf(form,tol,r,n,d)

tol = 1.0e-11 r = 3 + 1/(7 + 1/(16 + 1/(-294 + 1/(3 + 1/(-4))))) = 1146408/364913

#### Latex

If we're careful with the curly braces, Latex and MathJax can typeset that last continued fraction and produce a graphic for this post.

{{1}\over{{7+ {{1}\over{{16+ {{1}\over{{-294+ {{1}\over{{3+ {{1}\over{-4}}}}}}}}}}}}}}

$${{1}\over{{7+ {{1}\over{{16+ {{1}\over{{-294+ {{1}\over{{3+ {{1}\over{-4}}}}}}}}}}}}}}$$

#### ratp

I have another function, `ratp`, that is `rat` with `round` replaced by `floor` so that it produces proper, but longer, continued fractions with positive terms. Let's see how it compares to `rat`. Since the MATLAB quantity `pi` that approximates $\pi$ is rational, we can set `tol` to zero without fear of an infinite loop.

form = ' tol = %6.1e \n r = %s \n p = %s \n\n'; for tol = 10.^[-2 -6 -8 -11 -inf] fprintf(form,tol,rat(pi,tol),ratp(pi,tol)) end

tol = 1.0e-02 r = 3 + 1/(7) p = 3 + 1/(7) tol = 1.0e-06 r = 3 + 1/(7 + 1/(16)) p = 3 + 1/(7 + 1/(15 + 1/(1))) tol = 1.0e-08 r = 3 + 1/(7 + 1/(16 + 1/(-294))) p = 3 + 1/(7 + 1/(15 + 1/(1 + 1/(292)))) tol = 1.0e-11 r = 3 + 1/(7 + 1/(16 + 1/(-294 + 1/(3 + 1/(-4))))) p = 3 + 1/(7 + 1/(15 + 1/(1 + 1/(292 + 1/(1 + 1/(1 + 1/(1 + 1/(2)))))))) tol = 0.0e+00 r = 3 + 1/(7 + 1/(16 + 1/(-294 + 1/(3 + 1/(-4 + 1/(5 + 1/(-15))))))) p = 3 + 1/(7 + 1/(15 + 1/(1 + 1/(292 + 1/(1 + 1/(1 + 1/(1 + 1/(2 + 1/(1 + 1/(3 + 1/(1 + 1/(14))))))))))))

#### phi

The world's most irrational number, as measured by the length of continued fraction approximations, is $\phi$, the golden ratio. (Proof: solve $x = 1 + 1/x$.) The coefficients in the proper continued fraction are all ones, thereby achieving the slowest possible convergence. For the same accuracy, `ratp(phi,tol)` requires about twice as many terms as `rat(phi,tol)`.

```
format short
phi = (1 + sqrt(5))/2
tol = 1.0e-4
r = rat(phi)
p = ratp(phi)
```

phi = 1.6180 tol = 1.0000e-04 r = '2 + 1/(-3 + 1/(3 + 1/(-3 + 1/(3 + 1/(-3 + 1/(3 + 1/(-3)))))))' p = '1 + 1/(1 + 1/(1 + 1/(1 + 1/(1 + 1/(1 + 1/(1 + 1/(1 + 1/(1 + 1/(1 + 1/(1 + 1/(1 + 1/(1 + 1/(1 + 1/(1))))))))))))))'

tol = 0; lengthr = length(rat(phi,tol)) lengthp = length(ratp(phi,tol))

lengthr = 154 lengthp = 297

#### DIY

(For my international readers, DIY stands for "Do It Yourself".) Create your own `ratp` by making a one-word edit in a copy of `toolbox\matlab\specfun\rat.m`. Compare the two functions on `sqrt(2)`. Warning: all is not what appears at first glance.

Get the MATLAB code

Published with MATLAB® R2018b

**29**views (last 30 days) | |

**Category:**- Fun,
- History,
- Numerical Analysis,
- Precision

### Comments

To leave a comment, please click here to sign in to your MathWorks Account or create a new one.

## Recent Comments