Discrepancy with np.roots.

[Boomer91]

import numpy as np
from fqs import single_quartic

p = [0.1, 0.0, 0.0, 0.0, -1.0]
sol = single_quartic(*p)
sol_numpy = np.roots(p)
print(np.sort(sol))
print(np.sort(sol_numpy))

Returns:

[-3.16227766+0.j          0.        -3.16227766j -0.        +3.16227766j       3.16227766-0.j ]
[-1.77827941e+00+0.j          1.14491749e-16-1.77827941j      1.14491749e-16+1.77827941j      1.77827941e+00+0.j ]

Which is quite different. I suspect this is an edge case because of the many zeros. It can be worked out exaclty by hand,
but I was wondering why the code does not produce the correct answer.

I also noticed a few if s == 0:, would it not be more robust to use np.isclose(s, 0.0)?

Thanks,
Cyriel

[mzivic7]

Solver used in this project has few cases where it generates large errors:
Cardano's method for solving depressed cubic has large round-off error
Ferrari's method for solving depressed quartic becomes computationally unstable when largest real solution to Ferrari's resolvent cubic (z0 in the code) approaches zero.

These problems are fixed by modifying Cardano's and Ferrari's algorithms as described here:
https://quarticequations.com

Here are some results:

Ferraris mod: [0.+1.77827941j  0.-1.77827941j  1.77827941+0.j  -1.77827941-0.j]
np.polyroots: [-3.33066907e-16+1.77827941j  -3.33066907e-16-1.77827941j  1.77827941e+00+0.j  -1.77827941e+00+0.j]
fqs:          [-0.+3.16227766j  0.-3.16227766j  3.16227766-0.j  -3.16227766+0.j]
wolframalpha: [ 0.+1.77828j  0.+-1.77828j  1.77828+0.j  -1.77828+0.j]