When using the `roots`

function from the *Polynomials* package, the values found for multiple roots can be quite inaccurate. Just as an example, here is the polynomial with zeros at 1 and -1, each of multiplicity 4.

```
julia> using Pynomials
julia> p = Poly([1, 0, -4, 0, 6, 0, -4, 0, 1])
Poly(1 - 4*x^2 + 6*x^4 - 4*x^6 + x^8)
julia> Polynomials.roots(p)
8-element Array{Complex{Float64},1}:
-1.000063312800363 + 6.33123618252134e-5im
-1.000063312800363 - 6.33123618252134e-5im
-0.9999366871996364 + 6.33132485906425e-5im
-0.9999366871996364 - 6.33132485906425e-5im
1.0000782706569455 + 0.0im
0.9999999979609331 + 7.826862289848257e-5im
0.9999999979609331 - 7.826862289848257e-5im
0.9999217334211853 + 0.0im
```

I understand that this inaccuracy is introduced by the eigenvalue function that finds all these roots in one go. Each of these roots can be refined by applying, e.g., Newtonâs or Laguerreâs method.

I thought, package *PolynomialRoots* would provide such a feature but, at least in the way I use it, returns a similarly inaccurate solution.

```
julia> using PolynomialRoots
julia> p = complex([1.0, 0, -4, 0, 6, 0, -4, 0, 1]);
julia> PolynomialRoots.roots(p)
8-element Array{Complex{Float64},1}:
-1.0000703710890786 + 6.890751979805258e-9im
0.9999999994387893 + 8.914437556930486e-5im
-1.0000000085541791 - 7.037275037868812e-5im
0.9999999926924392 - 8.91443743617057e-5im
1.0000891483097958 - 3.3737785193170825e-9im
-0.999929625585618 - 6.892300583852365e-9im
0.9999108595589759 + 3.3725709200610236e-9im
-0.9999999947711241 + 7.037275192729228e-5im
```

Is there a way â or a Julia function â to âpolishâ these roots all in once and without utilizing multiple-precision numbers?