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

[Dan Sebald]

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

Reiterating Rik's comment, I guess that is pretty much the same as what I
pointed out.  The


nz = find (abs (r) > tol);


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:


        if (abs (ra (i_a) - rb (i_b)) < tol)


which tests if the roots of the different polynomials are the same.  That
could be changed to something like:


        if (abs (ra (i_a) - rb (i_b)) < tol*norm ([ra(i_a) rb(i_b) 1]))


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:


octave-cli:39> test polygcd
  ***** assert (polygcd (poly (1:8), poly (3:12)), poly (3:8), sqrt (eps))
!!!!! test failed
assert (polygcd (poly (1:8), poly (3:12)),poly (3:8),sqrt (eps)) expected
       1     -33     445   -3135   12154  -24552   20160
but got
 Columns 1 through 6:

   1.0000e+00  -3.3000e+01   4.4500e+02  -3.1350e+03   1.2154e+04 
-2.4552e+04

 Column 7:

   2.0160e+04
maximum absolute error 9.27736e-05 exceeds tolerance 1.49012e-08


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

Additional Item Attachment:

File name: octave-polygcd_2014aug06_2.patch Size:1 KB


    _______________________________________________________

Reply to this item at:

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

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