Polynomial.Roots() does not return the correct results

[sfetzel]

var polynomialCoefficients = new double[]
{
   1.19469367e+21, 1.67924808e+16, 3.97850921e+10, 1.39471145e+00, 1.56417732e-07
};
var polynomial = new Polynomial(polynomialCoefficients);

The result of polynomial.Roots() will return the following roots:
[0]: {(-4458290,732664504, 504313068,0086204)}
[1]: {(-4458290,732664504, -504313068,0086204)}
[2]: {(0, 0)}
[3]: {(0, 0)}

Using octave:
roots([1.56417732e-07, 1.39471145e+00, 3.97850921e+10, 1.67924808e+16, 1.19469367e+21])
The expected roots are:
(-4247248.3831557, 5.04311305e+08),
(-4247248.3831557, -5.04311305e+08),
(-331498.88729568, 0),
(-90585.83564952, 0)

My current workaround is:

Evd<Complex> eigen = polynomial.EigenvalueMatrix().ToComplex().Evd(Symmetricity.Asymmetric);
var roots = eigen.EigenValues;

However, these roots are not 100% complex conjugated.

[Arlofin]

The implementation of Polynomial.Roots() is almost identical to the described workaround, with the difference that it stays in the real numbers by omitting ToComplex(). Algebraically, this is a valid algorithm as well, but for the example polynomial it seems to be numerically less stable.
I do not have much experience with the topic, but intuitively, using complex numbers is never a bad idea for solving algebraic equations. So my suggestion would be to insert a call to ToComplex() in the implementation of Polynomial.Roots().
A further note on the example polynomial: This polynomial is ill-conditioned, hence numerical difficulties are to be expected. The complex solutions produced by the workaround might not be 100% conjugated, but they are nevertheless more precise than those computed with Octave (and the deviation from the conjugation symmetry [occurring at the 13th digit of the real part] is smaller than the deviation of the Octave roots from the true root [deviation at the 9th digit]). For very high precision computation of polynomial roots neither MathNet.Numerics nor Octave are the right tools as they are limited by concept. For those who are interested, I recommend to use Julia in combination with the PolynomialRoots package (or some other specialized packages) and BigFloat arithmetic.

[Arlofin]

Using the complex eigenvalue decomposition would, however, run into issue #595, and other issues as the unit tests hang up even when applying the fix in #944 .

[sfetzel]

I worked on an implementation of the algorithm explained in the paper "Fast and Backward Stable Computation of Roots of Polynomials" (see https://epubs.siam.org/doi/10.1137/140983434) for a university course and used C# for that with MathNet.Numerics. It should be more precise in theory. See the fork here: https://github.com/sfetzel/mathnet-numerics
Do you think it would make sense to incorporate this code into the library?

Nevertheless, I think the location of the AmvwSolver.cs I chose is incorrect. There are some things, which would need to be improved, for example: the computation of the Givens Rotations can be more precise.