1## Copyright (C) 2001 Paul Kienzle 2## 3## This program is free software; you can redistribute it and/or modify 4## it under the terms of the GNU General Public License as published by 5## the Free Software Foundation; either version 2 of the License, or 6## (at your option) any later version. 7## 8## This program is distributed in the hope that it will be useful, 9## but WITHOUT ANY WARRANTY; without even the implied warranty of 10## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11## GNU General Public License for more details. 12## 13## You should have received a copy of the GNU General Public License 14## along with this program; if not, write to the Free Software 15## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 16 17## -*- texinfo -*- 18## @deftypefn {Function File} {[@var{n}, @var{d}]} = rat (@var{x}, @var{tol}) 19## 20## Find a rational approximation to @var{x} within tolerance defined 21## by @var{tol} using a continued fraction expansion. E.g, 22## 23## @example 24## rat(pi) = 3 + 1/(7 + 1/16) = 355/113 25## rat(e) = 3 + 1/(-4 + 1/(2 + 1/(5 + 1/(-2 + 1/(-7))))) = 1457/536 26## @end example 27## 28## @end deftypefn 29## @seealso{rats} 30 31function [n,d] = rat(x,tol) 32 33 if (nargin != [1,2] || nargout != 2) 34 usage("[n,d] = rat(x,tol)"); 35 end 36 37 y = x(:); 38 39 ## replace inf with 0 while calculating ratios 40 y(isinf(y)) = 0; 41 42 ## default norm 43 if (nargin < 2) 44 tol = 1e-6 * norm(y,1); 45 end 46 47 ## First step in the approximation is the integer portion 48 n = round(y); # first element in the continued fraction 49 d = ones(size(y)); 50 frac = y-n; 51 lastn = ones(size(y)); 52 lastd = zeros(size(y)); 53 54 ## grab new factors until all continued fractions converge 55 while (1) 56 ## determine which fractions have not yet converged 57 idx = find (abs(y-n./d) >= tol); 58 if (isempty(idx)) break;end 59 60 ## grab the next step in the continued fraction 61 flip = 1./frac(idx); 62 step = round(flip); # next element in the continued fraction 63 frac(idx) = flip-step; 64 65 ## update the numerator/denominator 66 nextn = n; 67 nextd = d; 68 n(idx) = n(idx).*step + lastn(idx); 69 d(idx) = d(idx).*step + lastd(idx); 70 lastn = nextn; 71 lastd = nextd; 72 endwhile 73 74 ## move the minus sign to the top 75 n = n.*sign(d); 76 d = abs(d); 77 78 ## return the same shape as you receive 79 n = reshape(n, size(x)); 80 d = reshape(d, size(x)); 81 82 ## use 1/0 for Inf 83 n(isinf(x)) = sign(x(isinf(x))); 84 d(isinf(x)) = 0; 85 86 ## use 0/0 for NaN 87 n(isnan(x)) = 0; 88 d(isnan(x)) = 0; 89 90endfunction 91