Mon 10 Oct 2016 08:03:31 PM UTC, comment #11:
I deleted the code that was removing small coefficients. See http://hg.savannah.gnu.org/hgweb/octave/rev/dc872d5d74c5. This also may have fixed bug #34266. The fix was made on the development branch as we are very close to the 4.2.0 release.
|
Mon 10 Oct 2016 06:46:06 PM UTC, comment #10:
br has size 1x4, with br(1) around -9e-35 and br(2) around -4e-25.
|
Mon 10 Oct 2016 06:16:52 PM UTC, comment #9:
@Marco: Thanks for testing. Can you try another one from the BIST tests? This is from https://savannah.gnu.org/bugs/?34266.
If we eliminate all scrubbing code for small values then Octave returns a value for br which is 1x4 where the first two entries are nearly zero. If Matlab does the same thing then we can be reasonably certain they are not applying the filtering present in polyreduce.
There is still a possibility that they are replacing values less than eps(1) with 0. But first I want to see if Matlab returns a 1x4 vector for br.
|
Mon 10 Oct 2016 05:43:50 PM UTC, comment #8:
Uhm... matlab's residue does not accept 4 outputs. All the asserts (except that on e) are correct, length(br)=4 with br(1) exactly 0.
|
Mon 10 Oct 2016 04:41:04 PM UTC, comment #7:
@Marco: Could you try the following code in Matlab? It is from one of the BIST tests for residue. If Matlab is really not applying any filtering for small results than br will be a 4-element vector where the first entry is nearly 0.
|
Mon 10 Oct 2016 11:52:47 AM UTC, comment #6:
To comment #1
I think the code which tries to remove small coefficients should be removed, since too dangerous. I tried the following in matlab
b(2) is about 6.66e-16 and less than max([abs(a),abs(b)])*eps. But it was not converted to 0 by matlab.
|
Fri 07 Oct 2016 10:15:18 PM UTC, comment #5:
Opps. Never mind. I just checked. The no common roots.
|
Fri 07 Oct 2016 10:11:36 PM UTC, comment #4:
@RIk,
The polynomials "a" and "b" are each 3rd order. I haven't checked, but perhaps the r, p, & k you proposed produced a numerator which shares a root with the denominator?
|
Fri 07 Oct 2016 07:26:00 PM UTC, comment #3:
More thoughts, since residue is invertible it should be able to apply residue twice and arrive back at the original inputs. This works for more ordinary inputs, but not for the special case presented here. I made a script, tst_res.m, which is attached.
When running it in Octave I get,
As you can see, r2 and p2 don't even have the correct number of elements, and k2 is wildly off from 0.
(file #38683)
|
Fri 07 Oct 2016 04:49:42 PM UTC, comment #2:
As a second example,
As you can see, b2 is nearly equivalent to the starting variable b. However, if the code to trim polynomial coefficients was changed to act separately on the real part and on the imaginary part then these small numbers would be removed and there would be complete equality. Perhaps that is an addition to the code that is worth considering separately from the first issue reported.
|
Fri 07 Oct 2016 04:34:35 PM UTC, comment #1:
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.
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
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.
|
Fri 07 Oct 2016 12:12:48 PM UTC, original submission:
The Octave result
seems to be incorrect (the first element of den is missing). With MATLAB, I get what seems to be the correct result:
|