--- a/scripts/polynomial/mpoles.m	Thu Apr 08 16:04:13 2021 -0700
+++ b/scripts/polynomial/mpoles.m	Wed Apr 14 15:03:09 2021 +0200
@@ -29,7 +29,8 @@
 ## @deftypefnx {} {[@var{multp}, @var{idxp}] =} mpoles (@var{p}, @var{tol}, @var{reorder})
 ## Identify unique poles in @var{p} and their associated multiplicity.
 ##
-## The output is ordered from largest pole to smallest pole.
+## The output is ordered from pole with largest magnitude to smallest
+## magnitude.
 ##
 ## If the relative difference of two poles is less than @var{tol} then they are
 ## considered to be multiples.  The default value for @var{tol} is 0.001.
@@ -71,25 +72,19 @@
 
   Np = numel (p);
 
-  ## Force the poles to be a column vector.
+  ## force poles to be a column vector
 
   p = p(:);
 
-  ## Sort the poles according to their magnitidues, largest first.
-
   if (reorder)
-    ## Sort with smallest magnitude first.
-    [p, ordr] = sort (p);
-    ## Reverse order, largest maginitude first.
-    n = Np:-1:1;
-    p = p(n);
-    ordr = ordr(n);
+    ## sort with largest magnitude first
+    [~, ordr] = sort (abs (p), "descend");
+    p = p(ordr);
   else
-    ordr = 1:Np;
+    ordr = (1:Np).';
   endif
 
-  ## Find pole multiplicity by comparing the relative difference in the
-  ## poles.
+  ## find pole multiplicity by comparing relative difference of poles
 
   multp = zeros (Np, 1);
   indx = [];
@@ -124,3 +119,8 @@
 %!test
 %! [mp, n] = mpoles ([0 0], 0.01);
 %! assert (mp, [1; 2]);
+
+%!test
+%! [mp, n] = mpoles ([-1e4, -0.1, 0]);
+%! assert (mp, ones (3, 1));
+%! assert (n, [1; 2; 3]);

--- a/scripts/polynomial/residue.m	Thu Apr 08 16:04:13 2021 -0700
+++ b/scripts/polynomial/residue.m	Wed Apr 14 15:03:09 2021 +0200
@@ -360,11 +360,12 @@
 %! assert (b, [1, -5, 9, -3, 1], 1e-12);
 %! assert (a, [1, -5, 8, -4], 1e-12);
 %! [rr, pr, kr, er] = residue (b, a);
-%! [jnk, n] = mpoles (p);
-%! assert (rr, r(n), 1e-12);
-%! assert (pr, p(n), 1e-12);
+%! [~, m] = mpoles (rr);
+%! [~, n] = mpoles (r);
+%! assert (rr(m), r(n), 1e-12);
+%! assert (pr(m), p(n), 1e-12);
 %! assert (kr, k, 1e-12);
-%! assert (er, e(n), 1e-12);
+%! assert (er(m), e(n), 1e-12);
 
 %!test
 %! b = [1];
@@ -413,3 +414,12 @@
 %! k = 0;
 %! [p, q] = residue (r, pin, k);
 %! assert (p(4), 4.6828e+42, -1e-5);
+
+%!test <*60384>
+%! B = [1315.789473684211];
+%! A = [1, 1.100000536842105e+04, 1.703789473684211e+03, 0];
+%! poles1 = roots (A);
+%! [r, p, k, e] = residue (B, A);
+%! [B1, A1] = residue (r, p, k, e);
+%! assert (B, B1);
+%! assert (A, A1);