|
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 |
|
7016
|
7 ## the Free Software Foundation; either version 3 of the License, or (at |
|
|
8 ## your option) any later version. |
|
6788
|
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 |
|
7016
|
16 ## along with Octave; see the file COPYING. If not, see |
|
|
17 ## <http://www.gnu.org/licenses/>. |
|
6788
|
18 |
|
|
19 ## -*- texinfo -*- |
|
|
20 ## @deftypefn {Function File} {@var{s} =} rat (@var{x}, @var{tol}) |
|
|
21 ## @deftypefnx {Function File} {[@var{n}, @var{d}] =} rat (@var{x}, @var{tol}) |
|
|
22 ## |
|
|
23 ## Find a rational approximation to @var{x} within tolerance defined |
|
|
24 ## by @var{tol} using a continued fraction expansion. E.g, |
|
|
25 ## |
|
|
26 ## @example |
|
|
27 ## rat(pi) = 3 + 1/(7 + 1/16) = 355/113 |
|
|
28 ## rat(e) = 3 + 1/(-4 + 1/(2 + 1/(5 + 1/(-2 + 1/(-7))))) = 1457/536 |
|
|
29 ## @end example |
|
|
30 ## |
|
7001
|
31 ## Called with two arguments returns the numerator and denominator separately |
|
6788
|
32 ## as two matrices. |
|
|
33 ## @end deftypefn |
|
|
34 ## @seealso{rats} |
|
|
35 |
|
|
36 function [n,d] = rat(x,tol) |
|
|
37 |
|
|
38 if (nargin != [1,2] || nargout > 2) |
|
|
39 print_usage (); |
|
|
40 endif |
|
|
41 |
|
|
42 y = x(:); |
|
|
43 |
|
|
44 ## replace inf with 0 while calculating ratios |
|
|
45 y(isinf(y)) = 0; |
|
|
46 |
|
|
47 ## default norm |
|
|
48 if (nargin < 2) |
|
|
49 tol = 1e-6 * norm(y,1); |
|
|
50 endif |
|
|
51 |
|
|
52 ## First step in the approximation is the integer portion |
|
|
53 n = round(y); # first element in the continued fraction |
|
|
54 d = ones(size(y)); |
|
|
55 frac = y-n; |
|
|
56 lastn = ones(size(y)); |
|
|
57 lastd = zeros(size(y)); |
|
|
58 |
|
|
59 nd = ndims(y); |
|
6967
|
60 nsz = numel (y); |
|
6788
|
61 steps = zeros([nsz, 0]); |
|
|
62 |
|
|
63 ## grab new factors until all continued fractions converge |
|
|
64 while (1) |
|
|
65 ## determine which fractions have not yet converged |
|
|
66 idx = find (abs(y-n./d) >= tol); |
|
|
67 if (isempty(idx)) break; endif |
|
|
68 |
|
|
69 ## grab the next step in the continued fraction |
|
|
70 flip = 1./frac(idx); |
|
|
71 step = round(flip); # next element in the continued fraction |
|
|
72 |
|
|
73 if (nargout < 2) |
|
|
74 tsteps = NaN .* ones (nsz, 1); |
|
|
75 tsteps (idx) = step; |
|
|
76 steps = [steps, tsteps]; |
|
|
77 endif |
|
|
78 |
|
|
79 frac(idx) = flip-step; |
|
|
80 |
|
|
81 ## update the numerator/denominator |
|
|
82 nextn = n; |
|
|
83 nextd = d; |
|
|
84 n(idx) = n(idx).*step + lastn(idx); |
|
|
85 d(idx) = d(idx).*step + lastd(idx); |
|
|
86 lastn = nextn; |
|
|
87 lastd = nextd; |
|
|
88 endwhile |
|
|
89 |
|
|
90 if (nargout == 2) |
|
|
91 ## move the minus sign to the top |
|
|
92 n = n.*sign(d); |
|
|
93 d = abs(d); |
|
|
94 |
|
|
95 ## return the same shape as you receive |
|
|
96 n = reshape(n, size(x)); |
|
|
97 d = reshape(d, size(x)); |
|
|
98 |
|
|
99 ## use 1/0 for Inf |
|
|
100 n(isinf(x)) = sign(x(isinf(x))); |
|
|
101 d(isinf(x)) = 0; |
|
|
102 |
|
|
103 ## reshape the output |
|
|
104 n = reshape (n, size (x)); |
|
|
105 d = reshape (d, size (x)); |
|
|
106 else |
|
|
107 n = ""; |
|
|
108 nsteps = size(steps, 2); |
|
|
109 for i = 1: nsz |
|
|
110 s = [int2str(y(i))," "]; |
|
|
111 j = 1; |
|
|
112 |
|
|
113 while (true) |
|
|
114 step = steps(i, j++); |
|
|
115 if (isnan (step)) |
|
|
116 break; |
|
|
117 endif |
|
|
118 if (j > nsteps || isnan (steps(i, j))) |
|
|
119 if (step < 0) |
|
|
120 s = [s(1:end-1), " + 1/(", int2str(step), ")"]; |
|
|
121 else |
|
|
122 s = [s(1:end-1), " + 1/", int2str(step)]; |
|
|
123 endif |
|
|
124 break; |
|
|
125 else |
|
|
126 s = [s(1:end-1), " + 1/(", int2str(step), ")"]; |
|
|
127 endif |
|
|
128 endwhile |
|
|
129 s = [s, repmat(")", 1, j-2)]; |
|
|
130 n = cat (1, n, s); |
|
|
131 endfor |
|
|
132 endif |
|
|
133 |
|
|
134 endfunction |
