[Octave-bug-tracker] [bug #49291] residue function gives incorrect result

[Rik]

Follow-up Comment #1, bug #49291 (project octave):

This seems to be entirely due to a cleanup routine in residue.m.

I used the debugger to step through the code while executing your example, and
I have the Matlab-equivalent result until the end of the file.


pden =

   1.0000e+00   1.4481e+05   4.9290e+09   4.5817e+13   1.0008e+17


At that point, the m-file tries to remove elements which it thinks are likely
due to round-off errors in the calculation.  In this case, anything less than
the machine precision of the largest value.  The code starting at line 341 is


  ## Check for leading zeros and trim the polynomial coefficients.
  if (isa (r, "single") || isa (p, "single") || isa (k, "single"))
    small = max ([max(abs(pden)), max(abs(pnum)), 1]) * eps ("single");
  else
    small = max ([max(abs(pden)), max(abs(pnum)), 1]) * eps;
  endif

  pnum(abs (pnum) < small) = 0;
  pden(abs (pden) < small) = 0;

  pnum = polyreduce (pnum);
  pden = polyreduce (pden);

endfunction


The obvious question is whether this code is necessary.  If it is necessary,
should it apply both to the numerator (pnum) and the denominator (pden) or
just the numerator?  And if it is still necessary, and it should apply to both
pnum and pden, then how should the code be changed to recover the original
(and Matlab) result.

This code is really old and some form of this truncation has been in place
since 2007.  I'm adding the authors listed in the file to this bug report in
case they remember why this should be the way it is.

The Matlab documentation
(http://www.mathworks.com/help/matlab/ref/residue.html) says that the
denominator polynomial should be the size of 'r' or 'p' which in this case is
4.  So it is possible that this is also a bug/inconsistency in Matlab.


    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?49291>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/