Mercurial > octave
annotate scripts/general/rat.m @ 27919:1891570abac8
update Octave Project Developers copyright for the new year
In files that have the "Octave Project Developers" copyright notice,
update for 2020.
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Mon, 06 Jan 2020 22:29:51 -0500 |
| parents | b442ec6dda5c |
| children | bd51beb6205e |
| rev | line source |
|---|---|
|
27919
1891570abac8
update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents:
27918
diff
changeset
|
1 ## Copyright (C) 2001-2020 The Octave Project Developers |
|
27918
b442ec6dda5c
use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents:
26600
diff
changeset
|
2 ## |
|
b442ec6dda5c
use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents:
26600
diff
changeset
|
3 ## See the file COPYRIGHT.md in the top-level directory of this distribution |
|
b442ec6dda5c
use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents:
26600
diff
changeset
|
4 ## or <https://octave.org/COPYRIGHT.html/>. |
|
b442ec6dda5c
use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents:
26600
diff
changeset
|
5 ## |
| 6788 | 6 ## |
| 7 ## This file is part of Octave. | |
| 8 ## | |
|
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23573
diff
changeset
|
9 ## Octave is free software: you can redistribute it and/or modify it |
| 6788 | 10 ## under the terms of the GNU General Public License as published by |
|
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23573
diff
changeset
|
11 ## the Free Software Foundation, either version 3 of the License, or |
|
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22489
diff
changeset
|
12 ## (at your option) any later version. |
| 6788 | 13 ## |
| 14 ## Octave is distributed in the hope that it will be useful, but | |
| 15 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
|
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22489
diff
changeset
|
16 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
|
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22489
diff
changeset
|
17 ## GNU General Public License for more details. |
| 6788 | 18 ## |
| 19 ## You should have received a copy of the GNU General Public License | |
| 7016 | 20 ## along with Octave; see the file COPYING. If not, see |
|
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23573
diff
changeset
|
21 ## <https://www.gnu.org/licenses/>. |
| 6788 | 22 |
| 23 ## -*- texinfo -*- | |
| 26176 | 24 ## @deftypefn {} {@var{s} =} rat (@var{x}) |
| 25 ## @deftypefnx {} {@var{s} =} rat (@var{x}, @var{tol}) | |
| 26 ## @deftypefnx {} {[@var{n}, @var{d}] =} rat (@dots{}) | |
| 27 ## | |
| 28 ## Find a rational approximation of @var{x} to within the tolerance defined by | |
| 29 ## @var{tol}. | |
| 6788 | 30 ## |
| 26176 | 31 ## If unspecified, the default tolerance is @code{1e-6 * norm (@var{x}(:), 1)}. |
| 32 ## | |
| 33 ## When called with one output argument, return a string containing a | |
| 34 ## continued fraction expansion (multiple terms). | |
| 35 ## | |
| 36 ## When called with two output arguments, return numeric matrices for the | |
| 37 ## numerator and denominator of a fractional representation of @var{x} such | |
| 38 ## that @code{@var{x} = @var{n} ./ @var{d}}. | |
|
20158
7503499a252b
doc: Update docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
39 ## |
|
7503499a252b
doc: Update docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
40 ## For example: |
| 6788 | 41 ## |
| 42 ## @example | |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
9039
diff
changeset
|
43 ## @group |
|
26600
f6730533820e
doc: clean up doc example blocks for accuracy and consistent formatting
Mike Miller <mtmiller@octave.org>
parents:
26376
diff
changeset
|
44 ## s = rat (pi) |
| 26176 | 45 ## @result{} s = 3 + 1/(7 + 1/16) |
| 46 ## | |
|
26600
f6730533820e
doc: clean up doc example blocks for accuracy and consistent formatting
Mike Miller <mtmiller@octave.org>
parents:
26376
diff
changeset
|
47 ## [n, d] = rat (pi) |
|
f6730533820e
doc: clean up doc example blocks for accuracy and consistent formatting
Mike Miller <mtmiller@octave.org>
parents:
26376
diff
changeset
|
48 ## @result{} n = 355 |
|
f6730533820e
doc: clean up doc example blocks for accuracy and consistent formatting
Mike Miller <mtmiller@octave.org>
parents:
26376
diff
changeset
|
49 ## @result{} d = 113 |
| 26176 | 50 ## |
|
26600
f6730533820e
doc: clean up doc example blocks for accuracy and consistent formatting
Mike Miller <mtmiller@octave.org>
parents:
26376
diff
changeset
|
51 ## n / d - pi |
| 26176 | 52 ## @result{} 0.00000026676 |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
9039
diff
changeset
|
53 ## @end group |
| 6788 | 54 ## @end example |
| 55 ## | |
| 26176 | 56 ## Programming Note: With one output @code{rat} produces a string which is a |
| 57 ## continued fraction expansion. To produce a string which is a simple | |
| 58 ## fraction (one numerator, one denominator) use @code{rats}. | |
| 59 ## | |
| 60 ## @seealso{rats, format} | |
| 6788 | 61 ## @end deftypefn |
| 62 | |
| 26176 | 63 function [n, d] = rat (x, tol) |
| 6788 | 64 |
| 26176 | 65 if (nargin < 1 || nargin > 2) |
| 6788 | 66 print_usage (); |
| 67 endif | |
| 68 | |
| 26176 | 69 if (! isfloat (x)) |
| 70 error ("rat: X must be a single or double array"); | |
| 71 endif | |
| 72 | |
|
26197
95633ec174cf
rat.m: Add input validation to detect and fail on complex inputs (bug #55198)
Rik <rik@octave.org>
parents:
26176
diff
changeset
|
73 ## FIXME: This test should be removed when complex support is added. |
|
95633ec174cf
rat.m: Add input validation to detect and fail on complex inputs (bug #55198)
Rik <rik@octave.org>
parents:
26176
diff
changeset
|
74 ## See bug #55198. |
|
95633ec174cf
rat.m: Add input validation to detect and fail on complex inputs (bug #55198)
Rik <rik@octave.org>
parents:
26176
diff
changeset
|
75 if (iscomplex (x)) |
|
95633ec174cf
rat.m: Add input validation to detect and fail on complex inputs (bug #55198)
Rik <rik@octave.org>
parents:
26176
diff
changeset
|
76 error ("rat: X must be a real, not complex, array"); |
|
95633ec174cf
rat.m: Add input validation to detect and fail on complex inputs (bug #55198)
Rik <rik@octave.org>
parents:
26176
diff
changeset
|
77 endif |
|
95633ec174cf
rat.m: Add input validation to detect and fail on complex inputs (bug #55198)
Rik <rik@octave.org>
parents:
26176
diff
changeset
|
78 |
| 6788 | 79 y = x(:); |
| 80 | |
| 8506 | 81 ## Replace Inf with 0 while calculating ratios. |
| 26176 | 82 inf_idx = isinf (x); |
| 83 y(inf_idx(:)) = 0; | |
| 6788 | 84 |
| 26176 | 85 if (nargin == 1) |
| 86 ## default norm | |
| 87 tol = 1e-6 * norm (y, 1); | |
| 88 else | |
| 89 if (! (isscalar (tol) && isnumeric (tol) && tol > 0)) | |
| 90 error ("rat: TOL must be a numeric scalar > 0"); | |
| 91 endif | |
| 6788 | 92 endif |
| 93 | |
| 94 ## First step in the approximation is the integer portion | |
| 8506 | 95 |
| 96 ## First element in the continued fraction. | |
|
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14363
diff
changeset
|
97 n = round (y); |
|
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14363
diff
changeset
|
98 d = ones (size (y)); |
| 26176 | 99 frac = y - n; |
|
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14363
diff
changeset
|
100 lastn = ones (size (y)); |
|
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14363
diff
changeset
|
101 lastd = zeros (size (y)); |
| 6788 | 102 |
| 6967 | 103 nsz = numel (y); |
|
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14363
diff
changeset
|
104 steps = zeros ([nsz, 0]); |
| 6788 | 105 |
| 8506 | 106 ## Grab new factors until all continued fractions converge. |
| 6788 | 107 while (1) |
| 8506 | 108 ## Determine which fractions have not yet converged. |
| 26176 | 109 idx = find (y != 0 & abs (y - n./d) >= tol); |
|
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14363
diff
changeset
|
110 if (isempty (idx)) |
|
7881
f74669a09deb
rat.m: handle arrays and all-integer inputs
John W. Eaton <jwe@octave.org>
parents:
7031
diff
changeset
|
111 if (isempty (steps)) |
| 10549 | 112 steps = NaN (nsz, 1); |
|
7881
f74669a09deb
rat.m: handle arrays and all-integer inputs
John W. Eaton <jwe@octave.org>
parents:
7031
diff
changeset
|
113 endif |
|
f74669a09deb
rat.m: handle arrays and all-integer inputs
John W. Eaton <jwe@octave.org>
parents:
7031
diff
changeset
|
114 break; |
|
f74669a09deb
rat.m: handle arrays and all-integer inputs
John W. Eaton <jwe@octave.org>
parents:
7031
diff
changeset
|
115 endif |
| 6788 | 116 |
| 8506 | 117 ## Grab the next step in the continued fraction. |
| 26176 | 118 flip = 1 ./ frac(idx); |
| 8506 | 119 ## Next element in the continued fraction. |
|
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14363
diff
changeset
|
120 step = round (flip); |
| 6788 | 121 |
| 122 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
|
123 tsteps = NaN (nsz, 1); |
| 26176 | 124 tsteps(idx) = step; |
| 6788 | 125 steps = [steps, tsteps]; |
| 126 endif | |
| 127 | |
| 26176 | 128 frac(idx) = flip - step; |
| 6788 | 129 |
| 8506 | 130 ## Update the numerator/denominator. |
| 26176 | 131 savedn = n; |
| 132 savedd = d; | |
| 6788 | 133 n(idx) = n(idx).*step + lastn(idx); |
| 134 d(idx) = d(idx).*step + lastd(idx); | |
| 26176 | 135 lastn = savedn; |
| 136 lastd = savedd; | |
| 6788 | 137 endwhile |
| 138 | |
| 26176 | 139 if (nargout <= 1) |
| 140 ## string output | |
| 141 n = ""; | |
| 142 nsteps = columns (steps); | |
| 143 ## Loop over all values in array | |
| 144 for i = 1:nsz | |
| 145 | |
| 146 if (inf_idx(i)) | |
| 147 s = ifelse (x(i) > 0, "Inf", "-Inf"); | |
| 148 elseif (y(i) == 0) | |
| 149 s = "0"; | |
| 150 else | |
| 151 ## Create partial fraction expansion of one value | |
| 152 s = [int2str(y(i)), " "]; | |
| 153 j = 1; | |
| 154 | |
| 155 while (true) | |
| 156 step = steps(i, j++); | |
| 157 if (isnan (step)) | |
| 158 break; | |
| 159 endif | |
| 160 if (j > nsteps || isnan (steps(i, j))) | |
| 161 if (step < 0) | |
| 162 s = [s(1:end-1), " + 1/(", int2str(step), ")"]; | |
| 163 else | |
| 164 s = [s(1:end-1), " + 1/", int2str(step)]; | |
| 165 endif | |
| 166 break; | |
| 167 else | |
| 168 s = [s(1:end-1), " + 1/(", int2str(step), ")"]; | |
| 169 endif | |
| 170 endwhile | |
| 171 s = [s, repmat(")", 1, j-2)]; | |
| 172 endif | |
| 173 | |
| 174 ## Append result to output | |
| 175 n_nc = columns (n); | |
| 176 s_nc = columns (s); | |
| 177 if (n_nc > s_nc) | |
| 178 s(:, s_nc+1:n_nc) = " "; | |
| 179 elseif (s_nc > n_nc && n_nc != 0) | |
| 180 n(:, n_nc+1:s_nc) = " "; | |
| 181 endif | |
| 182 n = cat (1, n, s); | |
| 183 endfor | |
| 184 else | |
| 185 ## numerator, denominator output | |
| 186 | |
| 187 ## Move the minus sign to the numerator. | |
|
20231
83792dd9bcc1
Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents:
20158
diff
changeset
|
188 n .*= sign (d); |
|
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14363
diff
changeset
|
189 d = abs (d); |
| 6788 | 190 |
| 26176 | 191 ## Return the same shape as the input. |
|
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14363
diff
changeset
|
192 n = reshape (n, size (x)); |
|
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14363
diff
changeset
|
193 d = reshape (d, size (x)); |
| 6788 | 194 |
| 8506 | 195 ## Use 1/0 for Inf. |
| 26176 | 196 n(inf_idx) = sign (x(inf_idx)); |
| 197 d(inf_idx) = 0; | |
| 6788 | 198 endif |
| 199 | |
| 200 endfunction | |
|
12795
9e7ebbaf69ff
codesprint: new tests for files in scripts/general directory
John W. Eaton <jwe@octave.org>
parents:
12185
diff
changeset
|
201 |
|
9e7ebbaf69ff
codesprint: new tests for files in scripts/general directory
John W. Eaton <jwe@octave.org>
parents:
12185
diff
changeset
|
202 |
| 26176 | 203 %!assert (rat (pi), "3 + 1/(7 + 1/16)") |
| 204 %!assert (rat (pi, 1e-2), "3 + 1/7") | |
| 205 ## Test exceptional values | |
| 206 %!assert (rat (0), "0") | |
| 207 %!assert (rat (Inf), "Inf") | |
| 208 %!assert (rat (-Inf), "-Inf") | |
| 209 | |
|
12795
9e7ebbaf69ff
codesprint: new tests for files in scripts/general directory
John W. Eaton <jwe@octave.org>
parents:
12185
diff
changeset
|
210 %!test |
|
9e7ebbaf69ff
codesprint: new tests for files in scripts/general directory
John W. Eaton <jwe@octave.org>
parents:
12185
diff
changeset
|
211 %! [n, d] = rat ([0.5, 0.3, 1/3]); |
|
9e7ebbaf69ff
codesprint: new tests for files in scripts/general directory
John W. Eaton <jwe@octave.org>
parents:
12185
diff
changeset
|
212 %! assert (n, [1, 3, 1]); |
|
9e7ebbaf69ff
codesprint: new tests for files in scripts/general directory
John W. Eaton <jwe@octave.org>
parents:
12185
diff
changeset
|
213 %! assert (d, [2, 10, 3]); |
| 26176 | 214 ## Test exceptional values |
| 215 %!test | |
| 216 %! [n, d] = rat ([Inf, 0, -Inf]); | |
| 217 %! assert (n, [1, 0, -1]); | |
| 218 %! assert (d, [0, 1, 0]); | |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
219 |
|
23573
1b4f4ec53b4a
use new script to tag fixed bugs in tests
John W. Eaton <jwe@octave.org>
parents:
23572
diff
changeset
|
220 %!assert <*43374> (eval (rat (0.75)), [0.75]) |
|
19283
e1713e281ac5
rat.m: Remove unnecessary blank line at start of output (bug #43374).
Markus Bergholz <markuman+octave@gmail.com>
parents:
17744
diff
changeset
|
221 |
| 26176 | 222 ## Test input validation |
|
21581
6fab85c1538f
maint: Follow Octave conventions for use of semicolon in BIST tests.
Rik <rik@octave.org>
parents:
20852
diff
changeset
|
223 %!error rat () |
|
6fab85c1538f
maint: Follow Octave conventions for use of semicolon in BIST tests.
Rik <rik@octave.org>
parents:
20852
diff
changeset
|
224 %!error rat (1, 2, 3) |
| 26176 | 225 %!error <X must be a single or double array> rat (int8 (3)) |
|
26197
95633ec174cf
rat.m: Add input validation to detect and fail on complex inputs (bug #55198)
Rik <rik@octave.org>
parents:
26176
diff
changeset
|
226 %!error <X must be a real, not complex, array> rat (1+1i) |
| 26176 | 227 %!error <TOL must be a numeric scalar> rat (1, "a") |
| 228 %!error <TOL must be a numeric scalar> rat (1, [1 2]) | |
| 229 %!error <TOL must be a numeric scalar . 0> rat (1, -1) |