|
6788
|
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); |
|
6967
|
61 nsz = numel (y); |
|
6788
|
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 |
