function [N,D] = rat(X,tol)
%RAT	Rational approximation.
%	[N,D] = RAT(X,tol) returns two integer matrices so that N./D
%	is close to X in the sense that abs(N./D - X) <= tol*abs(X).
%	The rational approximations are generated by truncating continued
%	fraction expansions.   tol = 1.e-6*norm(X(:),1) is the default.
%
%	RAT(X) or RAT(X,tol) displays the continued fraction representation.
%
%	The same algorithm, with the default tol, is used internally
%	by MATLAB for FORMAT RAT.
%
%	See also FORMAT, RATS.

%	Cleve Moler, 10-28-90, 12-27-91, 9-4-92.
%	Copyright (c) 1984-93 by The MathWorks, Inc.


% Approximate x by
%
%                              1
%         d1 + ----------------------------------
%                                 1
%              d2 + -----------------------------
%                                   1
%                   d3 + ------------------------
%                                      1
%                        d4 + -------------------
%                                        1
%                             d5 + --------------
%                                           1
%                                  d6 + ---------
%                                             1
%                                       d7 + ----
%                                             d8
%
% The d's are obtained by repeatedly picking off the integer part and 
% then taking the reciprocal of the fractional part.  The accuracy of
% the approximation increases exponentially with the number of terms
% and is worst when x = sqrt(2).  For x = sqrt(2), the error with k
% terms is about 2.68*(.173)^k, which is
%
%         1    4.6364e-01
%         2    8.0210e-02
%         3    1.3876e-02
%         4    2.4006e-03
%         5    4.1530e-04
%         6    7.1847e-05
%         7    1.2430e-05
%         8    2.1503e-06
%         9    3.7201e-07
%        10    6.4357e-08

if nargin < 2
    tol = 1.e-6*norm(X(finite(X)),1);
else
    if isstr(tol), error('Use ''format rat'' instead.'), end
end
if nargout > 0
    N = zeros(size(X));
    D = zeros(size(X));
end
for j = 1:prod(size(X))
    x = X(j);
    k = 0;
    C = [1 0; 0 1];    % [n(k) n(k-1); d(k) d(k-1)];
    while 1
        k = k + 1;
        d = round(x);
        x = x - d;
        C = [C*[d;1] C(:,1)];
        if nargout == 0
            d = int2str(d);
            if k == 1
                s = d;
            elseif k == 2
                s = [s ' + 1/(' d ')'];
            else
                p = min(find(s==')'));
                s = [s(1:p-1) ' + 1/(' d ')' s(p:p+k-3)];
            end
        end
        if (abs(C(1,1)/C(2,1) - X(j)) <= tol) | (x == 0) | ~finite(x)
           break
        end
        x = 1/x;
    end
    if nargout == 0
        disp(s)
    else
        N(j) = sign(C(2,1))*C(1,1);
        D(j) = abs(C(2,1));
    end
end
