[Octave-bug-tracker] [bug #42742] polygcd fails valid test

[Rik]

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

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:


bad = 0;
tol = 10*sqrt (eps);
for ii=1:1000
  p = (unique (randn (10, 1)) * 10).';
  p1 = p(3:end);
  p2 = p(1:end-2);
  obs = polygcd (poly (-p1), poly (-p2), tol);
  exp = poly (- intersect (p1, p2));
  if (! size_equal (obs, exp))
    bad++;
  endif
endfor 


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


[d, r] = deconv (b, a);
nz = find (abs (r) > tol);
if (isempty (nz))


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)
    _______________________________________________________

Additional Item Attachment:

File name: polybad.var                    Size:0 KB


    _______________________________________________________

Reply to this item at:

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

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