Let me tell you about MATLAB's controversial functionrat。

Contents

Classicrat

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

<> helprat

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.  = RAT(X) produces integer matrices A and B so that

A ./ B = RAT(X)

Some examples:

long T = hilb(6), X = inv(T)  = 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)

一个例子

Let's try that second example with modern MATLAB andformat rat。

formatshortd = 1:6 e = ones(size(d)), A = abs(d'*e - e'*d) formatshorteX = inv(A) formatratX

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 ofX你更喜欢哪个?第一,印有一个浮动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 forratandformat rattoday, but they're still there and, as I said, they're moderately controversial.

The controversial algorithm

The core ofratgenerates 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 befloor。That always leaves a positive fraction to reciprocate. But 40 years ago, I got clever and usedroundinstead offloor。That means it may take fewer terms to obtain a specified accuracy, but the continued fractions are, shall I say,unorthodox。让我们看看more examples.

pi

Let's generate increasingly accurate rational approximations to $\pi$. With two output arguments,ratunwinds the continued fraction to produce two integers whose ratio has the same value. The accuracy ofratis determined by an optional tolerance,老鼠(x, tol)。The default tolerance is1.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 isratwithroundreplaced byfloorso that it produces proper, but longer, continued fractions with positive terms. Let's see how it compares torat。Since the MATLAB quantitypithat approximates $\pi$ is rational, we can settolto zero without fear of an infinite loop.

form =' tol = %6.1e \n r = %s \n p = %s \n\n';fortol = 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 asrat(phi,tol)。

formatshortphi = (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 ownratpby making a one-word edit in a copy oftoolbox\matlab\specfun\rat.m。Compare the two functions onsqrt(2)。Warning: all is not what appears at first glance.

得到了MATLAB代码

Published with MATLAB® R2018b