home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / SPECFUN.DI$ / RAT.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  2.9 KB  |  95 lines
  1. function [N,D] = rat(X,tol)
  2. %RAT    Rational approximation.
  3. %    [N,D] = RAT(X,tol) returns two integer matrices so that N./D
  4. %    is close to X in the sense that abs(N./D - X) <= tol*abs(X).
  5. %    The rational approximations are generated by truncating continued
  6. %    fraction expansions.   tol = 1.e-6*norm(X(:),1) is the default.
  7. %
  8. %    RAT(X) or RAT(X,tol) displays the continued fraction representation.
  9. %
  10. %    The same algorithm, with the default tol, is used internally
  11. %    by MATLAB for FORMAT RAT.
  12. %
  13. %    See also FORMAT, RATS.
  14.  
  15. %    Cleve Moler, 10-28-90, 12-27-91, 9-4-92.
  16. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  17.  
  18.  
  19. % Approximate x by
  20. %
  21. %                              1
  22. %         d1 + ----------------------------------
  23. %                                 1
  24. %              d2 + -----------------------------
  25. %                                   1
  26. %                   d3 + ------------------------
  27. %                                      1
  28. %                        d4 + -------------------
  29. %                                        1
  30. %                             d5 + --------------
  31. %                                           1
  32. %                                  d6 + ---------
  33. %                                             1
  34. %                                       d7 + ----
  35. %                                             d8
  36. %
  37. % The d's are obtained by repeatedly picking off the integer part and 
  38. % then taking the reciprocal of the fractional part.  The accuracy of
  39. % the approximation increases exponentially with the number of terms
  40. % and is worst when x = sqrt(2).  For x = sqrt(2), the error with k
  41. % terms is about 2.68*(.173)^k, which is
  42. %
  43. %         1    4.6364e-01
  44. %         2    8.0210e-02
  45. %         3    1.3876e-02
  46. %         4    2.4006e-03
  47. %         5    4.1530e-04
  48. %         6    7.1847e-05
  49. %         7    1.2430e-05
  50. %         8    2.1503e-06
  51. %         9    3.7201e-07
  52. %        10    6.4357e-08
  53.  
  54. if nargin < 2
  55.     tol = 1.e-6*norm(X(finite(X)),1);
  56. else
  57.     if isstr(tol), error('Use ''format rat'' instead.'), end
  58. end
  59. if nargout > 0
  60.     N = zeros(size(X));
  61.     D = zeros(size(X));
  62. end
  63. for j = 1:prod(size(X))
  64.     x = X(j);
  65.     k = 0;
  66.     C = [1 0; 0 1];    % [n(k) n(k-1); d(k) d(k-1)];
  67.     while 1
  68.         k = k + 1;
  69.         d = round(x);
  70.         x = x - d;
  71.         C = [C*[d;1] C(:,1)];
  72.         if nargout == 0
  73.             d = int2str(d);
  74.             if k == 1
  75.                 s = d;
  76.             elseif k == 2
  77.                 s = [s ' + 1/(' d ')'];
  78.             else
  79.                 p = min(find(s==')'));
  80.                 s = [s(1:p-1) ' + 1/(' d ')' s(p:p+k-3)];
  81.             end
  82.         end
  83.         if (abs(C(1,1)/C(2,1) - X(j)) <= tol) | (x == 0) | ~finite(x)
  84.            break
  85.         end
  86.         x = 1/x;
  87.     end
  88.     if nargout == 0
  89.         disp(s)
  90.     else
  91.         N(j) = sign(C(2,1))*C(1,1);
  92.         D(j) = abs(C(2,1));
  93.     end
  94. end
  95.