Re: Anyone familiar with polygcd?

[Michael D. Godfrey]

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:end-2);
> %!   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:end-2);
   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.