mathmore: invalid roots for a quartic polynomial

[sbinet]

trying to find the roots of the following polynomial:

  a x^4 + b x^3 + c x^2 + d x + e = 0

a = 2.2206846808021337
b = 7.643281053997895
c = 8.831759446092846 
d = 3.880673545129404
e = 0.5724551380144077

I get the following roots:

>>> import ROOT
>>> list(ROOT.Math.Polynomial(
2.2206846808021337, 7.643281053997895,
8.831759446092846, 3.880673545129404,
0.5724551380144077).FindRoots())

[(-0.8609914206605772+0j),
 (-0.8609285681365977+0j),
 (-0.8599689425356434-3.142626022335893e-05j),
 (-0.8599689425356434+3.142626022335893e-05j)]

plotting that polynomial, we can easily notice the roots are incorrect:

using numpy.roots (which is using a companion matrix):

>>> import numpy as np
>>> np.roots([
2.2206846808021337, 7.643281053997895,
8.831759446092846, 3.880673545129404,
0.5724551380144077])

array([-1.34283005+0.00000000e+00j, -1.34282804+0.00000000e+00j,
       -0.37809989+1.00708006e-06j, -0.37809989-1.00708006e-06j])

we do get saner values.

[lmoneta]

Hi,
Thank you for reporting this bug !
It looks to be a problem in a routine developed by HEP people, but it never came to GSL, it is copied inside Mathmore,
see https://root.cern/doc/master/complex__quartic_8h_source.html

If I change the code in ROOT::Math::Polynomial::FindRoots to use the generic GSL function, gsl_poly_complex_solve, it
works fine. This uses an iterative procedure instead of the analytical solution. GSL does not provide a function to solve analytically quartic equations.
Probably we should change the code in Polynomial and forget about that quartic implementation until somebody fixes it.
I am not sure if it is used somewhere else (e.g. CLHEP)

Lorenzo

[sbinet]

hi Lorenzo,

I stumbled on this when implementing my own GPL-free polynomial roots-finding algorithm from first principles.
I don't know the specifics of ROOT/GSL code, but in my code go-hep/exp#15 (BSD-3 licensed) I first go through building the cubic resolvent and then follow the wikipedia litterature.
note that for this particular case, what my code gives isn't completely satisfactory, though.
reading a bit more litterature, it would seem going through solving 2 quadratics would lead to better results.
I haven't got to try it out, though.

hth,
-s

[lmoneta]

Hi,
I see, if you are expert and can find the bug in this function, https://root.cern/doc/master/complex__quartic_8h_source.html, it will be great. I am sure somebody else might use that code

[lmoneta]

Looking back at my emails, I found that the routine comes from the LHCb RICH code and it was proposed to GSL, but it was never merged there at the end.

[lmoneta]

The problem seems to be in the calculation of disc at line 151
https://root.cern/doc/master/complex__quartic_8h_source.html#l00151

If I uncomment these following lines

//       more numerical problems with this calculation of disc
//      double CR2 = 729 * rcub * rcub;
//      double CQ3 = 2916 * qcub * qcub * qcub;
//       disc = (CR2 - CQ3) / 2125764.0;

it seems top provide the correct result in your example.
It would be nice to have more test points to validate the code

[sbinet]

I had found this polynomial to test the branch of my code where R==0. (brute force, by using the equivalent of MINUIT to minimize R).

(I am explicitly not looking at GPL-licensed code not to "taint" my BSD-licensed code)