Mercurial > octave-nkf
annotate scripts/general/rat.m @ 10549:95c3e38098bf
Untabify .m scripts
| author | Rik <code@nomad.inbox5.com> |
|---|---|
| date | Fri, 23 Apr 2010 11:28:50 -0700 |
| parents | 952d4df5b686 |
| children | be55736a0783 |
| rev | line source |
|---|---|
| 8920 | 1 ## Copyright (C) 2001, 2007, 2008, 2009 Paul Kienzle |
| 6788 | 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 ## | |
|
9039
51dc9691f23f
Cleanup documentation files errors.texi, debug.texi, io.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
23 ## Find a rational approximation to @var{x} within the tolerance defined |
|
51dc9691f23f
Cleanup documentation files errors.texi, debug.texi, io.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
24 ## by @var{tol} using a continued fraction expansion. For example, |
| 6788 | 25 ## |
| 26 ## @example | |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
9039
diff
changeset
|
27 ## @group |
| 7031 | 28 ## rat(pi) = 3 + 1/(7 + 1/16) = 355/113 |
| 29 ## rat(e) = 3 + 1/(-4 + 1/(2 + 1/(5 + 1/(-2 + 1/(-7))))) | |
| 30 ## = 1457/536 | |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
9039
diff
changeset
|
31 ## @end group |
| 6788 | 32 ## @end example |
| 33 ## | |
| 7001 | 34 ## Called with two arguments returns the numerator and denominator separately |
| 6788 | 35 ## as two matrices. |
| 36 ## @end deftypefn | |
| 37 ## @seealso{rats} | |
| 38 | |
| 39 function [n,d] = rat(x,tol) | |
| 40 | |
| 41 if (nargin != [1,2] || nargout > 2) | |
| 42 print_usage (); | |
| 43 endif | |
| 44 | |
| 45 y = x(:); | |
| 46 | |
| 8506 | 47 ## Replace Inf with 0 while calculating ratios. |
| 6788 | 48 y(isinf(y)) = 0; |
| 49 | |
| 50 ## default norm | |
| 51 if (nargin < 2) | |
| 52 tol = 1e-6 * norm(y,1); | |
| 53 endif | |
| 54 | |
| 55 ## First step in the approximation is the integer portion | |
| 8506 | 56 |
| 57 ## First element in the continued fraction. | |
| 58 n = round(y); | |
| 6788 | 59 d = ones(size(y)); |
| 60 frac = y-n; | |
| 61 lastn = ones(size(y)); | |
| 62 lastd = zeros(size(y)); | |
| 63 | |
| 64 nd = ndims(y); | |
| 6967 | 65 nsz = numel (y); |
| 6788 | 66 steps = zeros([nsz, 0]); |
| 67 | |
| 8506 | 68 ## Grab new factors until all continued fractions converge. |
| 6788 | 69 while (1) |
| 8506 | 70 ## Determine which fractions have not yet converged. |
|
7881
f74669a09deb
rat.m: handle arrays and all-integer inputs
John W. Eaton <jwe@octave.org>
parents:
7031
diff
changeset
|
71 idx = find(abs (y-n./d) >= tol); |
|
f74669a09deb
rat.m: handle arrays and all-integer inputs
John W. Eaton <jwe@octave.org>
parents:
7031
diff
changeset
|
72 if (isempty(idx)) |
|
f74669a09deb
rat.m: handle arrays and all-integer inputs
John W. Eaton <jwe@octave.org>
parents:
7031
diff
changeset
|
73 if (isempty (steps)) |
| 10549 | 74 steps = NaN (nsz, 1); |
|
7881
f74669a09deb
rat.m: handle arrays and all-integer inputs
John W. Eaton <jwe@octave.org>
parents:
7031
diff
changeset
|
75 endif |
|
f74669a09deb
rat.m: handle arrays and all-integer inputs
John W. Eaton <jwe@octave.org>
parents:
7031
diff
changeset
|
76 break; |
|
f74669a09deb
rat.m: handle arrays and all-integer inputs
John W. Eaton <jwe@octave.org>
parents:
7031
diff
changeset
|
77 endif |
| 6788 | 78 |
| 8506 | 79 ## Grab the next step in the continued fraction. |
| 6788 | 80 flip = 1./frac(idx); |
| 8506 | 81 ## Next element in the continued fraction. |
| 82 step = round(flip); | |
| 6788 | 83 |
| 84 if (nargout < 2) | |
|
10540
952d4df5b686
Eliminate NaN*ones and Inf*ones constructs and just use Nan() and Inf()
Rik <code@nomad.inbox5.com>
parents:
9051
diff
changeset
|
85 tsteps = NaN (nsz, 1); |
| 6788 | 86 tsteps (idx) = step; |
| 87 steps = [steps, tsteps]; | |
| 88 endif | |
| 89 | |
| 90 frac(idx) = flip-step; | |
| 91 | |
| 8506 | 92 ## Update the numerator/denominator. |
| 6788 | 93 nextn = n; |
| 94 nextd = d; | |
| 95 n(idx) = n(idx).*step + lastn(idx); | |
| 96 d(idx) = d(idx).*step + lastd(idx); | |
| 97 lastn = nextn; | |
| 98 lastd = nextd; | |
| 99 endwhile | |
| 100 | |
| 101 if (nargout == 2) | |
| 8506 | 102 ## Move the minus sign to the top. |
| 6788 | 103 n = n.*sign(d); |
| 104 d = abs(d); | |
| 105 | |
| 8506 | 106 ## Return the same shape as you receive. |
| 6788 | 107 n = reshape(n, size(x)); |
| 108 d = reshape(d, size(x)); | |
| 109 | |
| 8506 | 110 ## Use 1/0 for Inf. |
| 6788 | 111 n(isinf(x)) = sign(x(isinf(x))); |
| 112 d(isinf(x)) = 0; | |
| 113 | |
| 8506 | 114 ## Reshape the output. |
| 6788 | 115 n = reshape (n, size (x)); |
| 116 d = reshape (d, size (x)); | |
| 117 else | |
| 118 n = ""; | |
| 119 nsteps = size(steps, 2); | |
| 120 for i = 1: nsz | |
| 121 s = [int2str(y(i))," "]; | |
| 122 j = 1; | |
| 123 | |
| 124 while (true) | |
| 10549 | 125 step = steps(i, j++); |
| 126 if (isnan (step)) | |
| 127 break; | |
| 128 endif | |
| 129 if (j > nsteps || isnan (steps(i, j))) | |
| 130 if (step < 0) | |
| 131 s = [s(1:end-1), " + 1/(", int2str(step), ")"]; | |
| 132 else | |
| 133 s = [s(1:end-1), " + 1/", int2str(step)]; | |
| 134 endif | |
| 135 break; | |
| 136 else | |
| 137 s = [s(1:end-1), " + 1/(", int2str(step), ")"]; | |
| 6788 | 138 endif |
| 139 endwhile | |
| 140 s = [s, repmat(")", 1, j-2)]; | |
|
7881
f74669a09deb
rat.m: handle arrays and all-integer inputs
John W. Eaton <jwe@octave.org>
parents:
7031
diff
changeset
|
141 n_nc = columns (n); |
|
f74669a09deb
rat.m: handle arrays and all-integer inputs
John W. Eaton <jwe@octave.org>
parents:
7031
diff
changeset
|
142 s_nc = columns (s); |
|
f74669a09deb
rat.m: handle arrays and all-integer inputs
John W. Eaton <jwe@octave.org>
parents:
7031
diff
changeset
|
143 if (n_nc > s_nc) |
| 10549 | 144 s(:,s_nc+1:n_nc) = " " |
|
7881
f74669a09deb
rat.m: handle arrays and all-integer inputs
John W. Eaton <jwe@octave.org>
parents:
7031
diff
changeset
|
145 elseif (s_nc > n_nc) |
| 10549 | 146 n(:,n_nc+1:s_nc) = " "; |
|
7881
f74669a09deb
rat.m: handle arrays and all-integer inputs
John W. Eaton <jwe@octave.org>
parents:
7031
diff
changeset
|
147 endif |
| 6788 | 148 n = cat (1, n, s); |
| 149 endfor | |
| 150 endif | |
| 151 | |
| 152 endfunction |