Re: Anyone familiar with polygcd?

[Rik]

On 08/05/2014 10:30 AM, John W. Eaton wrote:
> On 07/07/2014 01: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
>
> I started noticing this too.
>
> Should we try to eliminate random values from the test suite?  It seems
> bad to me to use random values in tests which could fail due to extreme
> values.
>
> OTOH, this test failure may indicate that the algorithm used for polygcd
> is not as good as it should be.
>
> jwe
>

Michael Godfrey filed a bug report with some more information
(https://savannah.gnu.org/bugs/?42742).

I did a little bit more digging as well and it seems to be a tolerance
issue in the algorithm itself.  polygcd can be called with three arguments
as polygcd (a,b, tol).  If no tolerance is specified then sqrt (eps) is
used.  For three runs of 1000 trials I see 17/22/10 failures.  If I supply
a tolerance of 10*sqrt(eps) then the number of failures goes down to
2/1/1.  Maybe there is a way to pick the tolerance based on some
characteristic of the input (like the norm) that would be better.

--Rik