Mercurial > hg > octave-nkf
annotate scripts/general/rat.m @ 9039:51dc9691f23f
Cleanup documentation files errors.texi, debug.texi, io.texi
Spellcheck
Stylecheck
| author | Rik <rdrider0-list@yahoo.com> |
|---|---|
| date | Sun, 22 Mar 2009 13:14:15 -0700 |
| parents | eb63fbe60fab |
| children | 1bf0ce0930be |
| 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 | |
| 7031 | 27 ## rat(pi) = 3 + 1/(7 + 1/16) = 355/113 |
| 28 ## rat(e) = 3 + 1/(-4 + 1/(2 + 1/(5 + 1/(-2 + 1/(-7))))) | |
| 29 ## = 1457/536 | |
| 6788 | 30 ## @end example |
| 31 ## | |
| 7001 | 32 ## Called with two arguments returns the numerator and denominator separately |
| 6788 | 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 | |
| 8506 | 45 ## Replace Inf with 0 while calculating ratios. |
| 6788 | 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 | |
| 8506 | 54 |
| 55 ## First element in the continued fraction. | |
| 56 n = round(y); | |
| 6788 | 57 d = ones(size(y)); |
| 58 frac = y-n; | |
| 59 lastn = ones(size(y)); | |
| 60 lastd = zeros(size(y)); | |
| 61 | |
| 62 nd = ndims(y); | |
| 6967 | 63 nsz = numel (y); |
| 6788 | 64 steps = zeros([nsz, 0]); |
| 65 | |
| 8506 | 66 ## Grab new factors until all continued fractions converge. |
| 6788 | 67 while (1) |
| 8506 | 68 ## 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
|
69 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
|
70 if (isempty(idx)) |
|
f74669a09deb
rat.m: handle arrays and all-integer inputs
John W. Eaton <jwe@octave.org>
parents:
7031
diff
changeset
|
71 if (isempty (steps)) |
|
f74669a09deb
rat.m: handle arrays and all-integer inputs
John W. Eaton <jwe@octave.org>
parents:
7031
diff
changeset
|
72 steps = NaN .* ones (nsz, 1); |
|
f74669a09deb
rat.m: handle arrays and all-integer inputs
John W. Eaton <jwe@octave.org>
parents:
7031
diff
changeset
|
73 endif |
|
f74669a09deb
rat.m: handle arrays and all-integer inputs
John W. Eaton <jwe@octave.org>
parents:
7031
diff
changeset
|
74 break; |
|
f74669a09deb
rat.m: handle arrays and all-integer inputs
John W. Eaton <jwe@octave.org>
parents:
7031
diff
changeset
|
75 endif |
| 6788 | 76 |
| 8506 | 77 ## Grab the next step in the continued fraction. |
| 6788 | 78 flip = 1./frac(idx); |
| 8506 | 79 ## Next element in the continued fraction. |
| 80 step = round(flip); | |
| 6788 | 81 |
| 82 if (nargout < 2) | |
| 83 tsteps = NaN .* ones (nsz, 1); | |
| 84 tsteps (idx) = step; | |
| 85 steps = [steps, tsteps]; | |
| 86 endif | |
| 87 | |
| 88 frac(idx) = flip-step; | |
| 89 | |
| 8506 | 90 ## Update the numerator/denominator. |
| 6788 | 91 nextn = n; |
| 92 nextd = d; | |
| 93 n(idx) = n(idx).*step + lastn(idx); | |
| 94 d(idx) = d(idx).*step + lastd(idx); | |
| 95 lastn = nextn; | |
| 96 lastd = nextd; | |
| 97 endwhile | |
| 98 | |
| 99 if (nargout == 2) | |
| 8506 | 100 ## Move the minus sign to the top. |
| 6788 | 101 n = n.*sign(d); |
| 102 d = abs(d); | |
| 103 | |
| 8506 | 104 ## Return the same shape as you receive. |
| 6788 | 105 n = reshape(n, size(x)); |
| 106 d = reshape(d, size(x)); | |
| 107 | |
| 8506 | 108 ## Use 1/0 for Inf. |
| 6788 | 109 n(isinf(x)) = sign(x(isinf(x))); |
| 110 d(isinf(x)) = 0; | |
| 111 | |
| 8506 | 112 ## Reshape the output. |
| 6788 | 113 n = reshape (n, size (x)); |
| 114 d = reshape (d, size (x)); | |
| 115 else | |
| 116 n = ""; | |
| 117 nsteps = size(steps, 2); | |
| 118 for i = 1: nsz | |
| 119 s = [int2str(y(i))," "]; | |
| 120 j = 1; | |
| 121 | |
| 122 while (true) | |
| 123 step = steps(i, j++); | |
| 124 if (isnan (step)) | |
| 125 break; | |
| 126 endif | |
| 127 if (j > nsteps || isnan (steps(i, j))) | |
| 128 if (step < 0) | |
| 129 s = [s(1:end-1), " + 1/(", int2str(step), ")"]; | |
| 130 else | |
| 131 s = [s(1:end-1), " + 1/", int2str(step)]; | |
| 132 endif | |
| 133 break; | |
| 134 else | |
| 135 s = [s(1:end-1), " + 1/(", int2str(step), ")"]; | |
| 136 endif | |
| 137 endwhile | |
| 138 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
|
139 n_nc = columns (n); |
|
f74669a09deb
rat.m: handle arrays and all-integer inputs
John W. Eaton <jwe@octave.org>
parents:
7031
diff
changeset
|
140 s_nc = columns (s); |
|
f74669a09deb
rat.m: handle arrays and all-integer inputs
John W. Eaton <jwe@octave.org>
parents:
7031
diff
changeset
|
141 if (n_nc > s_nc) |
|
f74669a09deb
rat.m: handle arrays and all-integer inputs
John W. Eaton <jwe@octave.org>
parents:
7031
diff
changeset
|
142 s(:,s_nc+1:n_nc) = " " |
|
f74669a09deb
rat.m: handle arrays and all-integer inputs
John W. Eaton <jwe@octave.org>
parents:
7031
diff
changeset
|
143 elseif (s_nc > n_nc) |
|
f74669a09deb
rat.m: handle arrays and all-integer inputs
John W. Eaton <jwe@octave.org>
parents:
7031
diff
changeset
|
144 n(:,n_nc+1:s_nc) = " "; |
|
f74669a09deb
rat.m: handle arrays and all-integer inputs
John W. Eaton <jwe@octave.org>
parents:
7031
diff
changeset
|
145 endif |
| 6788 | 146 n = cat (1, n, s); |
| 147 endfor | |
| 148 endif | |
| 149 | |
| 150 endfunction |