Wed 06 Aug 2014 10:52:29 PM UTC, comment #4:
Reiterating Rik's comment, I guess that is pretty much the same as what I pointed out. The
is the issue. The vector r can have a very wide range in terms of magnitude. Yes, perhaps applying norm in some way would do it.
An alternative test for the roots()based approach might be to incorporate the magnitude of the root. That is, in the patch I submitted previously is this test:
which tests if the roots of the different polynomials are the same. That could be changed to something like:
I put the '1' in there to prevent decreasing the tolerance to something potentially smaller than eps. This works, but again, the tests fail because of something like the following:
The solution is good, but because the poly() routine multiplies all those roots, the many small differences in the roots all accumulate in a multiplicative way. So, again, the tolerance of the results may need refinement, but also in some sense I wish the roots function produced some more accurate results too.
(file #31857)

Wed 06 Aug 2014 09:24:49 PM UTC, comment #3:
I've attached a patch that implements polygcd() in the way I described. That is, uses the roots() function and finds common roots. However, the first impression is that it too has the same sorts of tolerance problems.
I increased the polygcd tolerance (i.e., the input parameter) and then get the results expected, but the resulting polynomial coefficients then don't pass the tolerance test, off by some factor slightly greater than 10*sqrt(eps).
I'm not sure about this numerical stability comment in the documentation for plygcd(). It seems to me the bigger issue is that applying the tolerance in the way it is done is the issue. Taking a look at some of the tests in polygcd(), it seems to me that in one case a coefficient like 1 should be within 10sqrt(eps), but then also a coefficient like 20,000 should also be accurate with 10sqrt(eps). It's easily possible with polynomial orders of say 8 to have numbers in the range of 1 to 2e4. There's nothing that takes the magnitude of numbers into account. One would prefer a relative test, say 0.001%, rather than absolute in this case.
(file #31856)

Wed 06 Aug 2014 08:44:53 AM UTC, comment #2:
In the polygcd() routine is this line:
Are the zeros of r always assured to be at the front of the array? Perhaps so.
Anyway, is there another way of solving this instead of using deconvolution? The documentation states:
"This is equivalent to the polynomial found by multiplying together all the common roots."
It is pretty straightforward to compute the roots of both of the polynomials using the roots() function. Then if one were to put the roots in order there could be a fast method of finding common roots (as opposed to a factorial permutation comparison of all combinations of roots). Here's an illustration.
All that is needed is the short simple loop inside a loop that picks the common roots within some tolerance. (There is the intersect() command, but that might require equality, I'm not sure.)
The question is whether the proposed approach is any more robust in the case of polynomial coefficients being large.

Tue 05 Aug 2014 08:36:48 PM UTC, comment #1:
The problem can be partially alleviated by using a tolerance. The default is sqrt (eps) but I find that if I use 10 times that value the frequency of an error drops from an average of 16/1000 iterations to a mean of 1/1000.
Test code I used was:
One possibility is that the tolerance, which defaults to sqrt (eps), should be based on some characteristic of the input vectors to polygcd such as the norm.
The relevant lines in polygcd are
In cases where it fails the remainder vector r is very nearly zero, i.e., the remainder polynomial is very small. But there is one coefficient which just manages to exceed the tolerance so the algorithm continues for one more trip through the while loop and gets the wrong answer.
A sample failing set of variables p, p1, p2, obs, exp is attached to the report as polybad.var for those who want a test case.
(file #31842)

Sat 12 Jul 2014 01:04:56 PM UTC, original submission:
This was first reported on the Mainbtainters List, but it appears
to be a bug in polygcd.
On 07/07/2014 06:57 PM, Rik wrote:
> All,
>
> I'm getting occasional failures from the test code in polygcd. It seems to
> happen about 10% of the time according to the following code:
>
> for i = 1:100
> bm(i) = test ("polygcd");
> endfor
> sum (bm)
>
> The test that fails is
>
> %!test
> %! for ii=1:10
> %! p = (unique (randn (10, 1)) * 10).';
> %! p1 = p(3:end);
> %! p2 = p(1:end2);
> %! assert (polygcd (poly (p1), poly (p2)), poly ( intersect (p1, p2)),
> sqrt (eps));
> %! endfor
>
> I tried a few different random seeds to see if I could fix the value to
> something that would always pass, but no luck.
>
> The simplest thing is to make this an %!xtest which can occasionally fail.
> But if someone understands polygcd and could suggest a way to modify the
> test that would be preferable.
>
> The error I get is a dimensional mismatch:
>
> !!!!! test failed
> ASSERT errors for: assert (polygcd (poly (p1), poly (p2)),poly
> (intersect (p1, p2)),sqrt (eps))
>
> Location  Observed  Expected  Reason
> . O(1x1) E(1x7) Dimensions don't match
>
>
> Rik
Here is where this goes wrong:
If the line [d, r] = deconv (b, a); in polygcd returns a r vector whose first element
is zero, then the line a = r / r(1); in polygcd fails. At this point a ends up
as a = 1 and x = 1, and the assert fails as above. I do not right now see why a = 1
due to divide by zero, but this is the cause of the failure.
The following script will, for me, reliably produce the failure as described above:
for ii=1:1000
p = (unique (randn (10, 1)) * 10).';
p1 = p(3:end);
p2 = p(1:end2);
assert (polygcd (poly (p1), poly (p2)), poly ( intersect (p1, p2)), sqrt (eps));
endfor
Obviously, the line a = r / r(1); should not be executed if r(1) = 0. But, it is
not clear to me right now what should be done instead.
I hope this helps.
