Mercurial > octave-nkf

comparison scripts/general/rat.m @ 6788:c81a0f3f5a82

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-07-23 22:05:29 by dbateman]
author dbateman
date Mon, 23 Jul 2007 22:05:30 +0000
parents
children 65a28e9de0a5
comparison
equal deleted inserted replaced
6787:963a19576024 6788:c81a0f3f5a82
1 ## Copyright (C) 2001 Paul Kienzle
2 ##
3 ## This file is part of Octave.
4 ##
5 ## Octave is free software; you can redistribute it and/or modify it
6 ## under the terms of the GNU General Public License as published by
7 ## the Free Software Foundation; either version 2, or (at your option)
8 ## any later version.
9 ##
10 ## Octave is distributed in the hope that it will be useful, but
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 ## General Public License for more details.
14 ##
15 ## You should have received a copy of the GNU General Public License
16 ## along with Octave; see the file COPYING. If not, write to the Free
17 ## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
18 ## 02110-1301, USA.
19
20 ## -*- texinfo -*-
21 ## @deftypefn {Function File} {@var{s} =} rat (@var{x}, @var{tol})
22 ## @deftypefnx {Function File} {[@var{n}, @var{d}] =} rat (@var{x}, @var{tol})
23 ##
24 ## Find a rational approximation to @var{x} within tolerance defined
25 ## by @var{tol} using a continued fraction expansion. E.g,
26 ##
27 ## @example
28 ## rat(pi) = 3 + 1/(7 + 1/16) = 355/113
29 ## rat(e) = 3 + 1/(-4 + 1/(2 + 1/(5 + 1/(-2 + 1/(-7))))) = 1457/536
30 ## @end example
31 ##
32 ## Called with two arguments returns the numerator and deniminator seperately
33 ## as two matrices.
34 ## @end deftypefn
35 ## @seealso{rats}
36
37 function [n,d] = rat(x,tol)
38
39 if (nargin != [1,2] || nargout > 2)
40 print_usage ();
41 endif
42
43 y = x(:);
44
45 ## replace inf with 0 while calculating ratios
46 y(isinf(y)) = 0;
47
48 ## default norm
49 if (nargin < 2)
50 tol = 1e-6 * norm(y,1);
51 endif
52
53 ## First step in the approximation is the integer portion
54 n = round(y); # first element in the continued fraction
55 d = ones(size(y));
56 frac = y-n;
57 lastn = ones(size(y));
58 lastd = zeros(size(y));
59
60 nd = ndims(y);
61 nsz = prod (size (y));
62 steps = zeros([nsz, 0]);
63
64 ## grab new factors until all continued fractions converge
65 while (1)
66 ## determine which fractions have not yet converged
67 idx = find (abs(y-n./d) >= tol);
68 if (isempty(idx)) break; endif
69
70 ## grab the next step in the continued fraction
71 flip = 1./frac(idx);
72 step = round(flip); # next element in the continued fraction
73
74 if (nargout < 2)
75 tsteps = NaN .* ones (nsz, 1);
76 tsteps (idx) = step;
77 steps = [steps, tsteps];
78 endif
79
80 frac(idx) = flip-step;
81
82 ## update the numerator/denominator
83 nextn = n;
84 nextd = d;
85 n(idx) = n(idx).*step + lastn(idx);
86 d(idx) = d(idx).*step + lastd(idx);
87 lastn = nextn;
88 lastd = nextd;
89 endwhile
90
91 if (nargout == 2)
92 ## move the minus sign to the top
93 n = n.*sign(d);
94 d = abs(d);
95
96 ## return the same shape as you receive
97 n = reshape(n, size(x));
98 d = reshape(d, size(x));
99
100 ## use 1/0 for Inf
101 n(isinf(x)) = sign(x(isinf(x)));
102 d(isinf(x)) = 0;
103
104 ## reshape the output
105 n = reshape (n, size (x));
106 d = reshape (d, size (x));
107 else
108 n = "";
109 nsteps = size(steps, 2);
110 for i = 1: nsz
111 s = [int2str(y(i))," "];
112 j = 1;
113
114 while (true)
115 step = steps(i, j++);
116 if (isnan (step))
117 break;
118 endif
119 if (j > nsteps || isnan (steps(i, j)))
120 if (step < 0)
121 s = [s(1:end-1), " + 1/(", int2str(step), ")"];
122 else
123 s = [s(1:end-1), " + 1/", int2str(step)];
124 endif
125 break;
126 else
127 s = [s(1:end-1), " + 1/(", int2str(step), ")"];
128 endif
129 endwhile
130 s = [s, repmat(")", 1, j-2)];
131 n = cat (1, n, s);
132 endfor
133 endif
134
135 endfunction