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