Mercurial > octave

changeset 33481:4d037eddd28f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
rat.m: Fix various corner cases and make output compatible with eval (bug #55198) Various corner cases with Inf or NaN inputs were failing with the previous code. E.g. `[n, d] = rat (complex (0, inf))` would give `n = Nan - Nani` and `d = 0`. Passing complex inputs with NaN would give an error about tolerance instead, though real inputs with NaN were handled properly. This patch fixes those various corner cases. Further, the string output for the above corner cases were failing to give back the original input when passed to `eval`, due to e.g. `0 + (Inf)*i` giving back `NaN + Infi`. (This is bug #31974, closed as Won't Fix). The output format of `rat()` for complex inputs has therefore been changed to return `complex (a, b)` instead of `(a) + (b) * i`. This new output string can be safely passed to `eval` in all cases. * rat.m: Fix various complex corner cases with Inf or NaN in them. Cosmetic change to return `complex (a, b)` instead of `(a) + (b) * i`, for compatibility with `eval`. Add and update BISTs for the above changes. Update documentation.
author Arun Giridhar <arungiridhar@gmail.com>
date Mon, 29 Apr 2024 10:08:44 -0400
parents 29282fcb0f2a
children 787937817a64 31ced69982f1
files scripts/general/rat.m
diffstat 1 files changed, 82 insertions(+), 25 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/general/rat.m	Sun Apr 28 21:27:05 2024 -0700
+++ b/scripts/general/rat.m	Mon Apr 29 10:08:44 2024 -0400
@@ -61,7 +61,7 @@
 ## @example
 ## @group
 ## s = rat (0.5 + i * pi)
-## @result{} s = (1 + 1/(-2)) + (3 + 1/(7 + 1/16)) * i
+## @result{} s = complex (1 + 1/(-2), 3 + 1/(7 + 1/16))
 ##
 ## [n, d] = rat (0.5 + i * pi)
 ## @result{} n =  113 + 710i
@@ -72,9 +72,14 @@
 ## @end group
 ## @end example
 ##
-## Programming Note: With one output @code{rat} produces a string which is a
-## continued fraction expansion.  To produce a string which is a simple
-## fraction (one numerator, one denominator) use @code{rats}.
+## Programming Notes:
+##
+## 1. With one output @code{rat} produces a string which is a continued
+## fraction expansion.  To produce a string which is a simple fraction
+## (one numerator, one denominator) use @code{rats}.
+##
+## 2. The string output produced by @code{rat} can be passed to @code{eval}
+## to get back the original input up to the tolerance used.
 ##
 ## @seealso{rats, format}
 ## @end deftypefn
@@ -89,6 +94,43 @@
     error ("rat: X must be a single or double array");
   endif
 
+  if (iscomplex (x))
+    if (nargout == 2)  # return numerator and denominator
+      if (nargin == 2)
+        [nr, dr] = rat (real (x), tol);
+        [ni, di] = rat (imag (x), tol);
+      else
+        [nr, dr] = rat (real (x));
+        [ni, di] = rat (imag (x));
+      endif
+
+      ## For inputs with inf, the output is set to 1/0 or -1/0.
+      ## Override that to +inf/1 or -inf/1.
+      ii = (dr == 0 & nr > 0); dr(ii) = 1; nr(ii) = +inf;
+      ii = (dr == 0 & nr < 0); dr(ii) = 1; nr(ii) = -inf;
+      ii = (di == 0 & ni > 0); di(ii) = 1; ni(ii) = +inf;
+      ii = (di == 0 & ni < 0); di(ii) = 1; ni(ii) = -inf;
+
+      d = lcm (dr, di);  # now this should always be nonzero
+      n = complex (nr .* (d ./ dr), ni .* (d ./ di));
+    elseif (nargout <= 1)  # string output
+      if (nargin == 2)
+        realstr = rat (real (x), tol);
+        imagstr = rat (imag (x), tol);
+      else
+        realstr = rat (real (x));
+        imagstr = rat (imag (x));
+      endif
+
+      nn = rows (realstr);
+      start  = repmat ("complex (", nn, 1);
+      mid    = repmat (", ",        nn, 1);
+      finish = repmat (")",         nn, 1);
+      n = [start, realstr, mid, imagstr, finish];
+    endif
+    return
+  endif
+
   y = x(:);
 
   ## Replace Inf with 0 while calculating ratios.
@@ -104,20 +146,6 @@
     endif
   endif
 
-  if (iscomplex (x))
-    if (nargout == 2)  # return numerator and denominator
-      [nr, dr] = rat (real (x), tol);
-      [ni, di] = rat (imag (x), tol);
-      d = lcm (dr, di);
-      n = nr .* (d ./ dr) + ni .* (d ./ di) * i;
-    elseif (nargout <= 1)  # string output
-      realstr = rat (real (x), tol);
-      imagstr = rat (imag (x), tol);
-      n = [repmat("(", rows (realstr), 1), realstr, repmat(") + (", rows (realstr), 1), imagstr, repmat(") * i", rows (imagstr), 1)];
-    endif
-    return
-  endif
-
   ## First step in the approximation is the integer portion
 
   ## First element in the continued fraction.
@@ -246,7 +274,7 @@
 
 ## Test complex scalar input
 %!test <*55198>
-%! assert (rat (complex (0.5, pi)), "(1 + 1/(-2)) + (3 + 1/(7 + 1/16)) * i");
+%! assert (rat (complex (0.5, pi)), "complex (1 + 1/(-2), 3 + 1/(7 + 1/16))");
 %! [n, d] = rat (complex (0.5, pi));
 %! assert (n, 113 + 710*i);
 %! assert (d, 226);
@@ -260,21 +288,50 @@
 %! assert (d, [887313, 49387, 49387, 887313]);
 %! assert (all (abs (n ./ d - x) <= 2e-6));
 %! str = rat (x);
-%! assert (str(1, :), "(0 + 1/(3 + 1/(4 + 1/(4 + 1/(4 + 1/4))))) + (1 + 1/(-20 + 1/(-2 + 1/(-3 + 1/(-6))))) * i");
-%! assert (str(2, :), "(-1 + 1/(5 + 1/(4 + 1/(4 + 1/4)))       ) + (1 + 1/(-2 + 1/(-2 + 1/(-3 + 1/8)))    ) * i");
-%! assert (str(3, :), "(-1 + 1/(5 + 1/(4 + 1/(4 + 1/4)))       ) + (-1 + 1/(2 + 1/(2 + 1/(3 + 1/(-8))))   ) * i");
-%! assert (str(4, :), "(0 + 1/(3 + 1/(4 + 1/(4 + 1/(4 + 1/4))))) + (-1 + 1/(20 + 1/(2 + 1/(3 + 1/6)))     ) * i");
+%! assert (str(1, :), "complex (0 + 1/(3 + 1/(4 + 1/(4 + 1/(4 + 1/4)))), 1 + 1/(-20 + 1/(-2 + 1/(-3 + 1/(-6)))))");
+%! assert (str(2, :), "complex (-1 + 1/(5 + 1/(4 + 1/(4 + 1/4)))       , 1 + 1/(-2 + 1/(-2 + 1/(-3 + 1/8)))    )");
+%! assert (str(3, :), "complex (-1 + 1/(5 + 1/(4 + 1/(4 + 1/4)))       , -1 + 1/(2 + 1/(2 + 1/(3 + 1/(-8))))   )");
+%! assert (str(4, :), "complex (0 + 1/(3 + 1/(4 + 1/(4 + 1/(4 + 1/4)))), -1 + 1/(20 + 1/(2 + 1/(3 + 1/6)))     )");
 
 ## Test complex exceptional inputs
 %!test <*55198>
-%! assert (rat (complex (inf, 0)), "(Inf) + (0) * i");
-%! assert (rat (complex (0, inf)), "(0) + (Inf) * i");
+%! assert (rat (complex (inf, 0)), "complex (Inf, 0)");
+%! assert (rat (complex (0, inf)), "complex (0, Inf)");
+%! assert (rat (complex (-inf, 0)), "complex (-Inf, 0)");
+%! assert (rat (complex (0, -inf)), "complex (0, -Inf)");
+%! assert (rat (complex (nan, 0)), "complex (NaN , 0)");
+%! assert (rat (complex (0, nan)), "complex (0, NaN )");
+
+%!test <*55198>
+%! [n, d] = rat (complex (inf, 0));
+%! assert (n, complex (inf, 0));
+%! assert (d, 1);
+%! [n, d] = rat (complex (0, inf));
+%! assert (n, complex (0, inf));
+%! assert (d, 1);
+%! [n, d] = rat (complex (-inf, 0));
+%! assert (n, complex (-inf, 0));
+%! assert (d, 1);
+%! [n, d] = rat (complex (0, -inf));
+%! assert (n, complex (0, -inf));
+%! assert (d, 1);
+%! [n, d] = rat (complex (nan, 0));
+%! assert (n, complex (nan, 0));
+%! assert (d, 1);
+%! [n, d] = rat (complex (0, nan));
+%! assert (n, complex (0, nan));
+%! assert (d, 1);
 
 ## Test eval with complex inputs
 %!test <*55198>
 %! x = complex (0.5, pi);
 %! assert (eval (rat (x)), x, 1e-6 * norm (x, 1))
 
+## Test eval with inf*i
+%!test <*55198>
+%! x = complex (0, inf);
+%! assert (eval (rat (x)), x, 1e-6 * norm (x, 1))
+
 %!assert <*43374> (eval (rat (0.75)), [0.75])
 
 ## Test input validation